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COMPUTER PROGRAM FOR DESIGN OF TWO-DIMENSIONAL SUPERSONIC 
TURBINE ROTOR BLADES WITH BOUNDARY-LAYER CORRECTION 
by Louis J. Goldman and Vincent J. Scullin 
Lewis Research Center 

SUMMARY 

A FORTRAN IV computer program for the design of two-dimensional supersonic 
rotor blade sections corrected for boundary-layer displacement thickness is presented. 
The ideal rotor is designed by the method of characteristics to produce vortex flow 
within the blade passage. The boundary-layer parameters are calculated by Cohen and 
Reshotko's method for laminar flow and Sasman and Cresci's method for turbulent flow. 
The program input consists essentially of the blade surface Mach number distribution 
and total flow conditions. The primary output is the corrected blade profile and the 
boundary-layer parameters. 


INTRODUCTION 

Methods for the design of supersonic turbines for possible use in turbopump and 
open-cycle auxiliary power systems, where high-energy fluids are used and high pres- 
sure ratios are available, have recently become of interest. The method of character- 
istics as applied to the two-dimensional isentropic flow of a perfect gas can be used for 
the design of the supersonic blading. Computer programs for the isentropic design of 
two-dimensional supersonic nozzles and rotor blade sections have been described in 
references 1 and 2, respectively. A computer program for the design of two- 
dimensional supersonic nozzles with boundary -layer correction is described in refer- 
ence 3. 

This report presents a computer program for the design of supersonic rotor blades 
where losses are accounted for by correcting the rotor profile for boundary -layer dis- 
placement thickness. The ideal rotor blade profile is obtained by using the computer 
program described in reference 2. Boundary-layer parameters are calculated by using 
the computer program described in reference 4. The final rotor blade profile is then 



obtained by adding the displacement thicknesses to the ideal blade coordinates. The 
program described herein is essentially a modified combination of the two programs 
described in references 2 and 4. 

The boundary-layer parameters (displacement and momentum thicknesses) are also 
used to calculate the conditions downstream of the rotor after the flow has mixed to a 
uniform state. The procedure described in reference 5 is used for this purpose. 

This report presents a description of the input and output and a complete FOR- 
TRAN IV listing of the program. A brief description of the computer program and 
method of design is also given. An example of the program output is included to indicate 
the use of the program and the results obtainable. 


METHOD OF ANALYSIS 

The design of two-dimensional supersonic rotor blades that are corrected for 
boundary-layer displacement thickness is described herein. The ideal rotor passage is 
designed by the method of characteristics as applied to the isentropic flow of a perfect 
gas. Boundary-layer parameters (displacement and momentum thicknesses) are then 
calculated for the ideal blades. The final rotor blade profile is obtained by adding the 
displacement thicknesses to the ideal rotor coordinates. 


Rotor Blade Description and Design 

The design of the ideal blade passage is based on establishing vortex flow within the 
passage by a procedure analogous to that given in reference 2. A typical passage 
is shown in figure 1. The passage consists essentially of three major parts: (1) inlet 
transition arcs, (2) circular arcs, and (3) outlet transition arcs. The inlet transition 
arcs (upper and lower) are required to convert the assumed uniform parallel flow at the 
passage inlet into vortex flow. The concentric circular arcs turn and maintain the vor- 
tex flow condition. The outlet transition arcs reconvert the vortex flow into uniform 
parallel flow at the passage exit. Straight-line segments parallel to the inlet and outlet 
flow direction complete the ideal passage. 

As seen from figure 1, the ideal passage is designed so that the outlet spacing is 
less than the inlet spacing. This is necessary if the corrected passage (i. e. , the pas- 
sage corrected by the boundary -layer displacement thicknesses) is to have equal outlet 
and inlet spacing (see fig. 2). For an ideal passage designed for impulse conditions 
(equal inlet and outlet Mach numbers), this is accomplished by having less circular 
turning for the outlet portion of the passage. That is, the outlet circular arcs JK and 
DE (fig. 1) are less than the corresponding inlet circular arcs IJ and CD. An iterative 


2 



Characteristic 


Vortex 

transition 




Figure 2. - Design of supersonic rotor blade section. 








procedure performed by the computer program determines the condition for which the 
corrected passage has equal inlet and outlet spacings. 


Boundary-Layer Calculations 

The boundary-layer parameters (displacement and momentum thicknesses) are 
calculated for the ideal rotor passage by using the computer program described in 4. 

The program uses Cohen and Reshotko's method (ref. 6) for the calculation of laminar 
boundary layers and Sasman and Cresci's method (ref. 7) for turbulent boundary layers. 
Curvature effects are not considered in these calculations. 

In the laminar regime, a single ordinary differential equation (the momentum in- 
tegral equation) is solved numerically. The results of this method (as explained in 
ref. 6) have to be extended for flows in highly favorable pressure gradients as might 
occur over some portions of the rotor blade. For turbulent flow, coupled first-order 
ordinary differential equations (the momentum and moment -of -momentum integral equa- 
tions) are solved using Runge-Kutta techniques. 

The displacement thicknesses obtained from the program are then added to the ideal 
rotor blade coordinates to obtain the corrected blade profile. As discussed previously, 
an iterative procedure is performed by the computer program to determine the condi- 
tions for which the corrected passage has equal inlet and outlet spacings. It should be 
noted that it is not always possible to obtain a blade design by this procedure because of 
boundary-layer separation or large boundary-layer growth. 

A number of options related to the boundary-layer calculations are also available to 
the user. Because of the adverse pressure gradients that exist within the rotor passage, 
the computer program will generally predict laminar separation at these locations. The 
program allows for the reattachment of the flow and continuation of the calculations for 
turbulent flow, if this is desired. The user can also, if he wishes, force transition to 
turbulent flow at any point of the calculations, including the inlet. 


Aftermixing Conditions 

The displacement and momentum thickness at the rotor exit (station 1, fig. 2) can 
be used to calculate the aftermixing conditions downstream of the rotor assuming that 
the flow mixes to a uniform state. Application of the continuity, momentum, and energy 
equations between stations 1 and 2 (fig. 2) results in the determination of aftermixing 
Mach number, flow angle, pressure ratio, and kinetic energy loss. The calculation 
procedure has been described in reference 5. 
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Subsonic and supersonic aftermixing axial Mach number solutions are possible for 
this loss model when the free-stream axial Mach number at the blade exit (before mix- 
ing) is supersonic. The subsonic solution corresponds to mixing plus oblique shock 
losses, whereas the supersonic solution corresponds to shockless mixing. A more de- 
tailed discussion of the different solutions can be found in reference 8. 


DESCRIPTION OF INPUT 


A description of the input for the FORTRAN IV computer program is given in this 
section. The input consists primarily of rotor inlet and outlet Mach number, upper- 
and lower- surface Mach number, inlet flow angle, specific -heat ratio, and total flow 
conditions. Either U. S. customary units or the International System of units may be 
used for input. This option is controlled by the input variable KEM. The program gas 
properties are set up for air. For gases other than air the changes required to the pro- 
gram are described in appendix A. 

The input format is shown in table I. The input variables are 


BETAN inlet flow angle, 0., deg 

VIN inlet Prandtl-Meyer angle, deg 

VLOW lower-surface Prandtl-Meyer angle, y^, deg 

VUP upper-surface Prandtl-Meyer angle, y , deg 

VOUT outlet Prandtl-Meyer angle, y , deg 

BETAT initial estimate of outlet flow angle, /S Q , deg (in the absence of a better es- 
timate use BETAN) 

DELV flow-turning increment (recommended value, 0. 1), Ay, deg 

(The flow-turning increment must be specified so that (y^ - y^)/Ay, 

(u Q - Vj)/Ay, (y u - yj)/Ay, and (y u - y Q )/Ay are all integers.) 

GAM specific -heat ratio, y 

NTURBU integer number of station on upper surface, if any, at which user wishes 
turbulent boundary layer to begin (If NTURBU is set equal to zero, the 
program will begin laminar boundary- layer calculations at the rotor inlet. 
Value of CTHETU must be specified; see section Instructions for Prepar- 
ing Input. NTURBU cannot be set equal to 1 because initial values of dis- 
placement and momentum thicknesses are not available to program. ) 


NTURBL same as NTURBU except for lower surface 
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XMAX 










NVP 


R 

PTZ 

TTZ 

XMAX 

TE 

CTHETU 

CTHETL 

KPRE 


KGRAD 

KSDE 

KLAM 

KMAIN 

KPROF 

KEM 


integer number of points desired in velocity profile of the boundary layer at 
each station (must have at least 1) 

gas constant, J/(kg)(K); (ft)(lbf)/(slug)(°R) 

2 2 

inlet relative total pressure, N/m ; lbf/ft 
inlet relative total temperature, K; °R 
blade chord, m; ft 

thickness of trailing edge (used only for aftermixing calculations), m; ft 

real variable indicating ratio of momentum thickness after reattachment to 
momentum thickness at laminar separation for upper surface 

same as CTHETU except for lower surface 

integer (0 or 1) indicating whether printing of output from PRECAL is de- 


sired: 

Output suppressed 0 

Output printed 1 


integer (0 or 1, see KPRE) indicating whether printing of surface gradients 
of velocity and Mach number is desired 

integer (0 or 1, see KPRE) indicating whether printing of solutions of lami- 
nar and turbulent differential equations is desired 

integer (0 or 1, see KPRE) indicating whether printing of laminar calcula- 
tions for location of instability and transition is desired 

integer (0 or 1, see KPRE) indicating whether printing of principal calculated 
boundary-layer parameters is desired 

integer (0 or 1, see KPRE) indicating whether printing of velocity profiles 
is desired 

integer (0 or 1) indicating which of two allowable sets of units are used in 


input: 

U. S. customary (pounds force, slugs, feet, seconds, 

degrees Rankine, and foot-pounds) 0 

International System (newtons, kilograms, meters, seconds, 

kelvin, and joules) 1 
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Instructions for Preparing Input 

Laminar separation and reattachment. - If NTURBU and NTURBL are set equal td 
zero, the program will begin laminar -boundary -layer calculations at the rotor inlet. 
Because of the adverse pressure gradients that exist on the upper and lower surfaces 
the program may predict separation on these surfaces. The program is set up so that 
reattachment of the flow is allowed to occur, and the calculations proceed for a turbulent 
boundary layer, if this is desired. Values for CTHETU and CTHETL (different from 
zero), which are the ratios of the momentum thickness after reattachment to the momen- 
tum thickness at laminar separation, must be specified as input. The values of CTHETU 
and CTHETL can be set equal to 1, which is essentially equivalent to assuming that 
transition to turbulent flow occurs at the point of imminent laminar separation. If 
CTHETU and CTHETL are set equal to zero, reattachment will not occur. 

Output. - Usually KPRE, KGRAD, KSDE, KLAM, and KPROF are set equal to zero. 
They do not give the main output of the boundary -layer section. However, if this addi- 
tional output is desired, these quantities are set equal to 1. A description of this output 
is given in appendix B. 


DESCRIPTION OF MAIN OUTPUT 

An example of the output from the program is shown in table n. The output is given 
in U. S. customary units, and each section has been numbered to correspond to the fol- 
lowing description: 

(1) Output 1 of the program is a listing of the input data used for the rotor design 
plus a listing of miscellaneous parameters including 

(a) Inlet, outlet, and surface critical velocity ratios M* and Mach numbers M. 

(b) Dimensionless blade spacing G*, chord C*, and solidity SIGMA. 

For each surface (lower and upper) items (2) to (5) are given for each iteration (if more 
than one iteration is required). 

(2) Output 2 gives the input needed by the boundary -layer program. The first two 
lines of this output are a listing of the input data except for NST (which is the total num- 
ber of station points) and UPMACH (which is the inlet relative Mach number). The re- 
maining output is the surface coordinates and velocity distribution: 

XOM dimensionless X-coordinate, X/C* 

YOM dimensionless Y-coordinate, Y/C* 

VVCR critical velocity ratio 

TWAL wall temperature, K; °R (wall temperature assumed equal to total 
temperature TTZ) 
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TABLE H. - EXAMPLE OF PROGRAM OUTPUT 


BETMIN) = 70.0000 DEC V(IN) = 39.1200 0E3 
DELTA / = 0.1000 3EG 

MISCELLANEOUS PARAMETERS 


M»( IN) = 

1.3257 

M(IM) = 2.4999 

K( OUT ) = 2.4999 

MiUOUT) = 

1 . 82 57 

M » ( LOW ) = 

1.5068 

M ( LDH ) * 1.7448 

M ( UP ) = 2.9672 

M* ( UP ) = 

L .9562 


INPJT BET A 1 0 ) = -61.1300 DFS 


GMIN) = 0.5944 



SIGMA (IN) * 

2. 7603 


Z • * 

1.6406 



G*I OUT I = 0.4210 



S IGMA(OUT) s 

3. 8966 


DESIGN PARAMETERS 

VIUP) = *9.1200 DEG V { OUT J = 39. 1200 DEG BETAIQUT) = -61.1 

Y(LDW) = 19.1200 DEG GAMMA = 1.4000 


•** UPPER SURFACE 


BOUNDARY 

LAYER - 

INPUT 


NST 

NVP 

GAM 

K 

64 

5 

1.400 

1714.70 


KPRE 

KGRAU 

KSDE 

KLAM 

0 

0 

0 

0 


XOM 

YOM 

VVCR 

TWAL 

-0 . 49 j 5 7 

-0. 43670 

1.62569 

1913.30 

-0.47711 

-0.38798 

1.82569 

1913. ?C 

-0.45365 

-0.33725 

1.82569 

1912.30 

-0.44013 

-0.28652 

1.82569 

1913.30 

-0.42172 

-0.23579 

1.02569 

1913. 30 

-0.4032b 

-0.1H506 

1. HZ569 

1913.30 

-0.38479 

-0.13433 

1.02569 

1913. 3C 

-0.36633 

-0.08361 

1.62569 

1913. 20 

-0.34787 

-0.03286 

1.R2569 

1915. 30 

-0.32943 

0.01785 

1. 62*69 

1913. 3G 

-0.31094 

0.06858 

1.82569 

1913.50 

-0.30783 

0.07688 

1.83951 

1913.30 

-0.30453 

0.08524 

1 . o5 31 fc 

1913. 30 

-0.30102 

0.09367 

1. ‘6663 

1913.30 

-0.29731 

0.10218 

l. 07993 

1913.30 

-0.29333 

0. 11080 

1.89305 

1913.30 

-0.28921 

0.11952 

1.9060! 

1913.30 

-0.28t7? 

0.12836 

1.91881 

1913. 3C 

-0.28019 

0.13734 

1.93143 

1913.30 

-0.27513 

0. 1 4b46 

1.94390 

1913.30 

-0.26985 

0.15580 

!. 95620 

1913. 30 

-0.25524 

0.17872 

1.95620 

191 3.30 

-0.23879 

0.20029 

1.95620 

1913. 30 

-0.22033 

0.22033 

1.95620 

1913.30 

-0.20029 

0.23870 

1.95620 

1913. 30 

-0.17872 

0.25524 

1.95620 

1913.30 

-0.15589 

0.26985 

1.95620 

1913.30 

-0.13 159 

0.28240 

1.95620 

1913.30 

-0.10657 

0.29280 

1.95620 

1913.30 

-0.08065 

0 .30098 

1.95620 

1913.30 

-0.05411 

0.30686 

1.95620 

1913. 30 

-0.0271b 

0.31041 

1.95620 

1913.30 

0.00009 

0.31160 

1.95620 

1913. 30 

0.0271b 

0.31041 

1. 95620 

1913. 30 

0.05411 

0. 30686 

1.95620 

1913. 20 

0.08065 

0.30098 

1.95620 

1913. 30 

0.10657 

0.29260 

1.95620 

1913. 30 

0.13169 

0.20240 

1./5620 

1913. 30 

0.15589 

0.^6985 

1. >5620 

1913, 30 

0.17872 

0.25524 

1. "*5620 

1913. 30 

0.20029 

0.^3870 

1.95620 

1913. 30 

0.22033 

0.22033 

I . 95620 

1913. JO 

0.23879 

0. 20029 

1. >5620 

1913.30 

0.24269 

0.19554 

1.95620 

191 3. 30 

0.24 >25 

0. 10715 

1. 74390 

1913.30 

0.25553 

0. 17889 

1.93143 

1913.30 

0.26159 

0. 1 7U74 

1.41081 

1913.30 

0.25732 

0.16260 

1.90601 

1913. 30 

0.27279 

0.15471 

1.59305 

1913. 30 

0.27800 

0. i46Bl 

l. 37993 

1913.30 

0.28293 

0. 1 3b96 

1. *^6o63 

1913.20 


PTZ 

676. oO 


KMA I N 
1 


Tl l 

1913.30 


KPR.OF 

0 


UPMACM 

2.4999 


XMAX 

0.041H3 


NTURB CTHET \ 

3 l.OOu 


KEM 

0 


>00 OEG 


r e 

O.UDO40 
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TABLE IT. - Continued. EXAMPLE OF PROGRAM OUTPUT 


0.2877V 

0.13117 

1.65316 

1913.30 

0.29229 

0.12342 

1.-83951 

1913.30 

0.29655 

0.11570 

1. *2569 

1913.30 

0.31712 

0.07b56 

1. d2 569 

1913.30 

0.33 763 

0.04142 

1.62569 

1913.30 

0. 35903 

0. 00428 

1.62569 

1913.30 

0.37855 

-0.03207 

1.82569 

1913.30 

0.3993V 

-0.07001 

1.H2569 

1913.30 

0.41952 

-0.10715 

1. J2569 

1913.30 

0.4V003 

-0.14430 

1.62569 

1913.30 

0 .46047 

-0.18144 

1.H2569 

1913.30 

0.48095 

-0.21858 

1. «2569 

1913.30 

0.50143 

-0.25572 

1 . 82569 

1913.30 


C PRELIMINARY CALCULATIONS 


PSZ 

= 

39.62063 

TSZ 

= 

030.41 14 

UZ 


3571.79850 

ASZ 

< 

1428.8039 

AT Z 

s 

2143.1355 

RHSZ 

= 

0. 2717091E-04 

RHTZ 

= 

0.2062952E-U3 

MUSZ 

= 

0.5372009E-06 

MUTZ 


0.9149355E-06 

NUSZ 

S 

0. 1977 L18E-01 

NUTZ 

= 

0.4433079E-02 

CP 

S 

6001.4498 

PR 

= 

0.67232 

TC 

= 

0.4892751E-02 

ARCL 


0.0737 


Station 

PRES 

UE 

MF 

PUPTZ 

NOVCR 

i 

39 . 6c063 

3571.79855 

2.499852 

0.053541 

1.325694 

2 

39 • 620b 3 

3571.79855 

2.499852 

0.058541 

1.825694 

3 

39.62063 

3571.79855 

2.499852 

0.053541 

1 .325694 

4 

39.6^063 

3571.79055 

2.499852 

0.058541 

1.325694 

5 

39.62063 

3571.79055 

2.499852 

0. 058541 

1.325694 

6 

39.62063 

3571.79655 

2.49985? 

0.058541 

1 .325694 

7 

39.6-2063 

3571.79055 

2.499052 

0.058541 

1.325694 

8 

39.62063 

3571.79855 

2.499852 

0.058541 

1 . 325694 

9 

39.62063 

3571.79855 

2.499852 

0.059541 

1. 325694 

10 

39.62063 

3571 .79056 

2.499852 

0. 059341 

i. 325694 

11 

39.62063 

3571.79854 

2.499052 

0.058541 

1.325694 

^ 12 

37.04862 

3590.83670 

2. 543043 

0.054741 

1.3 3951 4 

PRINCI PAL 

B3JN0ARY LAYER 

I NFOR u AT [ ON 





INSTABILITY ODES 

NOT OCCUR 





TRANSITION OCCURS AT STATION 

3 




SEPARATION DUES 

NUT OCCUR 





LAMINAR BOUNDARY 

L A Y £ - STA 

rroNs i to 2 




TURBULENT BOUNDARY LAYER - STATIONS 3 TJ 

64 





D I MENSIONLO 3L AI)t 

COORDINATES 


ANGLCS OF ROTATION 

r 


JNCORRECTTO 

CORRECT' 3 

translated 


station 

X 

Y 

YC3RR 

YT RAN 


1 

-9. 1,2073 

-C .01835 

-0.01345 

-0.03351 

70. OOUO 

2 

-0.0199'. 

-0.J1623 

-0.01763 

-0.03270 

70.0000 

3 

-0.01919 

-0 • (.1411 

-0.01523 

-0.03039 

70 .00 CO 

4 

-0.01641 

-0.01199 

-0.01325 

-0.02041 

7O.C0U0 

5 

-0. JL7r>u 

-0. 909B6 

-0.01125 

-0.OZ64O 

70.0000 

6 

-0.0 R. J / 

-0.0077'- 

-0. 009..4 

-0.0243 ) 

70 . 00 CO 

7 

-O.Ui&lO 

-0.00562 

-0.00721 

-0.02237 

70.0000 

8 

-o. oi 

-0.00360 

-0.00510 

-0. 02034 

70. CO.'O 

9 

-■'.014^5 

-0. 00135 

-0.0031? 

-0.01830 

70. 0000 

10 

-0.013?" 

0. C007*> 

-0.00111 

-3.01627 

70.00 00 

11 

-P.ul Jj 1 

0.C02U7 

0.00093 

-0.01417 

70. COCO 

12 

-0. jl 2." 8 

0.0032? 

0. 001 3 3 

-0. OL 381 

69. 00 JO 
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TABLE n. - Continued. EXAMPLE OF PROGRAM OUTPUT 


STATION 

X 

S 

OELSR 

THET 

DELTA 

F3RH 

FORM I 

1 

-3.020730 

0. 

-0. 

-0. 

-3. 

7. *563 

2.8586 

2 

-0.01*958 

0.002250 

0.000*79 

0.00006* 

0.000826 

7. *563 

2.6566 

3 

-0.019185 

0.00*516 

0.000386 

0. 000091 

0.03231* 

*.2**7 

1.4000 

* 

-3.018*13 

0.00677* 

0.000*33 

0.000099 

3.032331 

*. 3709 

1.4561 

5 

-3.0176*1 

0.0O9C33 

0.000*7* 

0.000106 

3.032397 

*.*55* 

1.4936 

6 

-3.016868 

0.011291 

0.000512 

0.000113 

0.002*90 

*. 5081 

1.5171 

7 

-0.016096 

0.0135*9 

0.0005*5 

0.000120 

0.032600 

*. 5393 

1.5309 

8 

-3.01532* 

0.015807 

0.000577 

0.000127 

3.032719 

*.5567 

1.5387 

9 

-3.01*551 

3.018865 

0.000607 

0.000133 

3.0328** 

*.5656 

1.5*26 

10 

-3.013779 

0.020323 

0. 000636 

0.000139 

0.032973 ■ 

*. 5691 

1.5*42 

11 

-3.013007 

0.022582 

0.0006*6 

0.0001*3 

0.033100 

*. 5296 

1.5266 

12 

-3.012876 

0.022952 

0.000669 

0.0001*6 

3.03332* 

*. 5825 

1.50*1 


STATION 

CF 

T AIM 

RTH 

□ TOY 

NUSS 

HTR AN 

LKN 

1 

3. 

0. 

-0. 

0. 

0. 

0. 

2.0^1 

2 

3.06568 

5.07*81 

3.0 

-261125.21 

3.38 

-1277. 6206 

2.021 

3 

3.00768 

1.331*4 

*.3 

-78*95.36 

2.69 

-384.0582 

0.608 

* 

3.00688 

1. 19193 

*. 7 

-70270.32 

3.51 

-3*3.8152 

0. 6J0 

5 

3.00636 

1.102*4 

5.0 

-6*994.51 

*• *6 

-318.0020 

0.608 

6 

3.00603 

1.0*491 

5.3 

-61632. 65 

5.28 

-301.406* 

0. 608 

7 

3.00581 

1.00701 

5.7 

-59368.07 

6.11 

-290.4732 

0 -juB 

8 

3.33566 

0.98103 

6.0 

-57836.59 

6.9* 

-282.9900 

0 . 5 «j8 

9 

3.00555 

0.96238 

6.3 

-56737.07 

7.7H 

-277.6003 

0. Ci 08 

10 

3.005*7 

0.9*835 

6.6 

-55909.82 

8.62 

-273.5528 

0.608 

11 

3.00559 

0.96822 

6.7 

-57081.50 

9.76 

-279.2856 

0.608 

12 

3.30572 

0.95860 

6.5 

-569*2.25 

9.77 

-278.60*3 

0.5 V6 


LOWER SURFACE **■ 


NTURB CTHtTA r E 

3 1.000 0.000*3 


XOM 

YOU 

WCR 

THAL 

-0.49557 

-0.076*2 

1.82569 

1913.30 

-0.46597 

-0.05078 

1.61169 

1913.30 

-0.47639 

-0.02653 

1. 79752 

1913.30 

-0.46685 

-0.00355 

1.78316 

1913.30 

-0.45735 

0.0162* 

1.76861 

1913.30 

-0.4*793 

0.03895 

1. 75388 

1913.30 

-0.438*9 

0.05865 

1. 73695 

1913.30 

-0.429H 

0.077*0 

1.72383 

1913.30 

-0.41979 

0.09530 

1.70850 

1913. 30 

-0.410*9 

0.11238 

1. ‘'9297 

1913.30 

-0.40123 

0.12871 

1.6772* 

1913. 30 

-0.39233 

0. 1**3* 

l.o6129 

1913. 30 

-0.382*1 

0. i. 5932 

1.6*511 

1913. 30 

-0.37363 

0.* 7369 

1.62871 

1913. 30 

-0.36*49 

0.107*9 

1.61208 

1 913. 30 

-0.35536 

0.20075 

1.59520 

1913.30 

-0.3*62* 

0.21350 

1. 57007 

1913.30 

-0.3371* 

0.22578 

1. 56069 

1913.30 

-0.32*35 

0.23761 

1. 5*303 

1913.30 

-0.31995 

0.2*902 

1.52509 

1913.30 

-0.30983 

0.26002 

1.50605 

1913. 30 

-0.2860* 

0.2 BoO* 

1. _»0685 

1913. 30 

-0.26032 

0.30983 

1. 50685 

1913.30 

-0.23232 

0.33136 

1.50685 

1913.30 

-0.20226 

0.35032 

1. 53685 

1913.30 

-0.17095 

0.36662 

1. 50685 

1913.30 

-0.13835 

0.38012 

l. 50605 

1913.30 

-0.10*73 

0. 39073 

1. 5068 5 

1913. 30 

-0.0702* 

0.39837 

1. 50685 

1913.30 

-0.03525 

0 • *0298 

1. 50685 

1913.30 

0.00033 

0. *0*5 1 

1.50685 

1913.30 

0.03525 

0. *0298 

1. 50685 

1913. 30 

0.0792* 

0. 39837 

1. 50605 

1913. 30 

0.10*73 

0.390/3 

1 . 506 d 5 

1913. 30 

0.13I35 

0.38012 

L. 50635 

1913.30 

0.17095 

0. 36662 

1. 33685 

1913.30 

0.20225 

0.35032 

1. '"•368 5 

1913.30 

0.23232 

0.33136 

1. L0685 

1913.30 

0.26032 

0.30938 

1. "0685 

1913. 3 r 

0.2b603 

0.33*69 

1. c 3685 

1913. 30 

0.27675 

0.29522 

1. 62 50 9 

1913. jO 

0.287*3 

0.2853^ 

1. “*303 

1913.30 

0.29833 

0.27507 

1 . * S069 

1913.30 

0.30913 

3.26*3* 

l. .730 / 

1913.30 

0.32c 15 

0 ■ £ 5 3 1 * 

1.-9520 

1913.30 

0.33122 

0.2414 j 

1.6120P 

1913.30 

0.3*233 

0.22923 

J. 62871 

1913.30 

0.35 iS5 

0.<U44 

l.r*51 1 

1913.30 

0.36635 

0.20306 

l.o6l29 

1913.30 

0.37653 

0.1390* 

1.6772* 

1913.30 


BOUNDARY LAYER - INPUT 


NS T 
60 


GAM 

1.403 


KPRE K3RAD 

0 0 


KSDE 

0 


PTZ 

676. 80 


KMAIN 

1 


rrz 

1913. 3D 


UPHACH 
2. *999 


KEH 

0 
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TABLE n. - Continued. 


EXAMPLE OF PROGRAM OUTPUT 


0.38825 0.17433 
0.4003b 0.13888 
0.4123V 0.1*264 
0.42*23 0.12556 
0.43453 0.10755 
0.44^37 0.08b55 
Q.4oi8? 0.06848 
0.47473 0.34724 
0.48793 0.02476 
0.50143 0.00091 


i. 69297 1913.30 

1.70850 1913.30 

1. 72363 1913.30 

1.73995 1923.30 

1. 75338 1913.30 

1.76861 1913.30 

1.73316 1913.30 

1. 79752 1913.30 

1.81169 1913.30 

1.82569 1913.30 


3S 


PRELIM] 

[NARY 

CALCULATIONS 

PSZ 

= 

39.62063 

TSZ 

= 

050. 4i 14 

UZ 

- 

3571.79858 

ASZ 

= 

1428.8039 

ATZ 

* 

2143.1355 

RHSZ 

. 

0.27 17091E-04 

RHTZ 

= 

0.2G62952E-03 

MUSZ 

* 

0.5372009E-06 

MUTZ 

= 

0.9149355E-06 

NUSZ 

= 

0. 1977118E-01 

NUTZ 

* 

0.4435079E-02 

CP 

* 

6001.4498 

PR 

= 

0.67232 

TC 

s 

0.4892751E-02 

ARCL 

= 

0.0592 


AT I ON 

PRES 

UE 

ME 

pop r z 

YOVCR 

1 

39.62063 

357L. 79855 

2. 499852 

u.038541 

1.325694 

2 

42.33229 

3544.40982 

2.457333 

0. 062548 

1.311694 

3 

45.18726 

3516.68137 

2.415483 

0. 066 766 

1.797521 

4 

48.19350 

3488.58292 

2.374237 

0.071208 

1.783159 

5 

51.35440 

3460.12845 

2. 333597 

0.075B78 

1.768615 

6 

54.67782 

3431.29196 

2.293511 

0.0B0789 

l. 753875 

7 

58.16675 

3402 .09058 

2.253987 

0.085944 

1.738949 

9 

61.82872 

3372.50253 

2.214980 

0.091354 

1.723825 

9 

65.66948 

3342.52209 

2. 176470 

0.097029 

1.708501 

10 

69.69449 

3312.14563 

2.138441 

0.102975 

1.S92975 

11 

73.91059 

3281.35953 

2.100862 

0.109206 

1.377239 

12 

78.32432 

3250.15289 

2.063711 

0.115727 

1.561287 


f PRINCIPAL BOUNDARY LAYER INFORMATION 


INSTABILITY DOES NOT OCCUR 
TRANSITION OCCURS AT STATION 3 
SEPARATION ODES NOT OCCUR 

LAMINAR BOUNDARY LAYER - STATIONS 1 TO 2 
TURBULENT BOUNDARY LAYER - STATIONS 3 TO 60 


DIMENSIONED BLADE COORDINATES ANCLES OF ROT A T I )N 10 




UNCORRECTED 

CORRECTED 

TRANSLATED 


STAT I ON 

X 

Y 

YCORR 

YTRAN 


1 

-0.02073 

-0.00320 

-0.30320 

-0.01 835 

70.00 00 

2 

-0.02033 

-0.00212 

-0.00121 

-3.01636 

6‘>. COi-O 

3 

-0.01993 

-o.ooiu 

-0.00043 

-3.01559 

68.00-iU 

4 

-0.01953 

-0.00015 

0.00356 

-0.01459 

67.0000 

5 

-0.01913 

0.00076 

0. 00150 

-3. 01366 

66.0000 

6 

-0 .01874 

0.00163 

0.00233 

-3.31278 

65 . OOv-D 

7 

-0. J1834 

0.00245 

0.03321 

-0.01194 

64. no'.u 

6 

-0.01795 

0.00324 

0.00400 

-0.01116 

63 . 00 CO 

9 

-0.01756 

0.00399 

0.00475 

-0.01040 

62. CO jO 

10 

-0.01717 

0.00470 

0.00547 

-3.00969 

61.00 JO 

11 

-0.UL678 

0.0053b 

0.03615 

-3.00900 

60.00 *0 

. 12 

-0.01640 

0.00604 

0.006*1 

-0.00635 

59.00-J0 
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TABLE H. - Concluded. EXAMPLE OF PROGRAM OUTPUT 


fSTATION 

X 

S 

OELSR 

THET 

DELTA 

FORM 

FORM I 

1 

-3.020730 

0 . 

-0. 

-0. 

-0. 

7.4563 

2.8536 

2 

-3,020328 

0.001.45 

0.003329 

0.000045 

3.000573 

7. 2629 

2.641 3 

3 

-D. 019927 

0.002236 

0.000254 

0.000063 

0.001525 

4.0559 

i . 4000 

4 

-3.019528 

3.G03277 

G. 000278 

0.000067 

0.031435 

4.1361 

1.4800 

5 

-3.31)131 

0.004271 

0.0C0299 

0. 000071 

3.001381 

4. 1941 

1 • 5509 

6 

-3.018733 

0.005223 

0. 00031' 

0.000075 

0.001340 

4. 2296 

1.6121 

7 

-3.018342 

0.006137 

0.000332 

0.000078 

0.001329 

4. 2462 

1.6 64 6 

8 

-3.017950 

0.007014 

0.000346 

0.000091 

3.031319 

4.2482 

1.7104 

9 

-3.017559 

0.007558 

0.000359 

0. 000085 

3.001313 

4. 2403 

1.7512 

10 

-3.017171 

0.008671 

0.000371 

0.000088 

0.001311 

4.2263 

i.7890 

11 

-3.016783 

3.009456 

0. 300393 

0.000091 

0.001312 

4. 2094 

1.8251 

L 12 

-3.016397 

0.010216 

0 . 000395 

0.000094 

0.001314 

4.1921 

1.8609 


STATION 

CF 

TAUW 

RTH 

DTDY 

NUSS 

HTR AN 

CRN 

1 

0. 

0 . 

-0. 

0 . 

0. 

0 . 

2.021 

2 

0.08363 

6.7 L- 206 

2.3 

-374773.63 

2.2b 

-1933.6740 

1 .867 

3 

3.0085? 

1.58470 

3. 3 

-91994.48 

1.51 

-450.0371 

O.Gjl 

4 

3.00749 

1.4.-355 

3. a 

-81970. 2D 

2. 14 

-401.0598 

0 . 6*t3 

? 

3.00664 

1.29963 

4.2 

-74224.09 

2.56 

-363. 1 600 

0.655 

6 

0.00599 

1.2G582 

4.7 

-68292.74 

2.93 

-334. 1 394 

0 . 

7 

0 . OD 5 A b 

l. 1 3438 

5.1 

-63699.37 

3.27 

-311.6676 

0.678 

8 

0.00508 

l .07854 

5.7 

-60037.27 

3.53 

-293.7474 

0. 690 

9 

3.00474 

1.03290 

6.2 

-369B5.8B 

3.88 

-278.8177 

0. 702 

10 

3.00445 

0 . 9 9340 

6 « B 

-54313.07 

4.15 

-265. 7403 

0.714 

11 

3.00419 

0.95733 

7.4 

-51849.90 

4.41 

-253.6890 

0. 7 t 6 

„ 12 

3.00395 

0.9 '221 

8.0 

-49473.16 

4.63 

-242.0599 

0. 7 .9 


= 1.207 XMX2 = 0.641 
*/tfCR2= 1.415 EBAR,=Q. 


AFTERMIXING PROPERTIES 
ROTOR WITH NO BOUNDARY LAYER CORRECTION 

XMFSl =2,499? SPACING =.310733 T F =0.00040 XM2 =1.5823 V///CR1 = 1.826 XMX1 

ALPHl = 61.133 AtPH2= 66.110 PT2/PT0= 0.2821 >T0/P2= 14.678 T2/lf3= 0.6663 

ROTOR WITH ttUU NO AK Y LAVES CORRECTION 

XMFSl =2.4??? SPACING =.015205 TE =0.00040 XM2 =1.8334 i/ZVCRl = 1. 325 XMXi = 1.237 XMX2 = 0.710 

ALPH1= 61.130 ALPH2= 67.221 PT2/PTC- 0.4543 PTO/P t = 13.310 T2/TT0= 0.5980 V/9CR2= 1.553 EBAR ?=0. 2 308 B E T A-N =0. 

SJPERSONIC SOLUTION 

XMFSl =2.4?)? SPACING =.315223 TE =0.03343 XM2 =2.4299 V/VCR1 = 1.426 XMX1 = 1.237 X MX 2 = 1.476 

ALPHl= 61.130 ALPH2= 52.603 PT2/PT0= 0.4847 PT0/P2= 31.602 T?/TTO= 0.45B5 Y/VCR?= 1.802 E6AR 2 = 0.1 3n64 E TA-N’O. 


ITERATION COMPLETED 

OETA(OJT) = -61.1300 DE(. 


■OUT) ♦ OcLirnTAL) 


G* (INI = 
G*( JUT) 


0.6 94h 
0.4210 


SIGMA (IN) = 2.7603 

SIGMA I OUT) = 3 r 966 


*01* UNIT05, EOF. 


REC= 03003 F I L = D0002 


6 ' 27 


7 o91 


8-. 34 
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(3) Output 3 lists the properties calculated for the upstream conditions: 


PSZ 

TSZ 

UZ 

ASZ (ATZ) 

RHSZ (RHTZ) 

MUSZ (MUTZ) 

NUSZ (NUTZ) 

CP 

PR 

TC 

ARCL 


2 2 

upstream static pressure, N/m ; lbf/ft 

upstream static temperature, K; °R 

upstream velocity, m/sec; ft/sec 

speed of sound based on upstream static (total) temperature, 
m/sec; ft/sec 

static (total) density based on upstream static (total) temper- 
ature, kg/m slug/ft^ 

dynamic viscosity based on upstream static (total) tempera- 
ture, (N)(sec)/m^; (lbf)(sec)/ft^ 

kinematic viscosity based on upstream static (total) temper- 
ature, m^/sec; ft^/sec 

specific heat at constant pressure, J/(kg)(K); (ft)(lbf)/ 
(slug)(°R) 

Prandtl number 

thermal conductivity, j/(m)(sec)(K); (ft)(lbf)/(ft)(sec)(°R) 

total distance along surface, m; ft 


The next part of this output gives the variables describing the flow along the surface: 


PRES 

static pressure, N/m^; lbf/ft^ 

UE 

free-stream velocity, m/sec; ft/sec 

ME 

free -stream Mach number 

POPTZ 

static -to-total pressure ratio 

VOVCR 

critical velocity ratio 


(4) Output 4 corresponds to KMAIN. It indicates the regions of laminar and turbu- 
lent boundary layers, and the stations at which instability, transition, and separation 
occur. It gives all the principal boundary-layer output parameters: 

X X-coordinate, m; ft 

S surface length, m; ft 

DELSR displacement thickness, m; ft 

THET momentum thickness, m; ft 
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DELTA boundary -layer thickness, m; ft 
FORM compressible form factor 
FORMI incompressible form factor 

The next part of the output gives the skin-friction and heat -transfer parameters: 


CF 

skin friction coefficient at wall 


TAUW 

shear stress at wall, N/m^; lbf/ft^ 


RTH 

momentum -thickness Reynolds number 


DTDY 

slope of temperature profile at wall, K/m; °R/ft 


NUSS 

local Nusselt number 


HTRAN 

2 2 
heat transfer per unit area, J/(sec)(m ); (ft)(lbf)/(sec)(ft ) 

CRN 

Reynolds analogy parameter 



(5) Output 5 gives the ideal and corrected surface coordinates. Translated surface 
coordinates are also provided so that coordinates for a blade profile may be obtained. 

(6) After the output for both surfaces is given, the final output, numbered 6, lists 
the aftermixing properties. Variable names ending in 0 refer to the upstream station, 
while those ending in 1 refer to the rotor exit (before mixing) and those ending in 2 refer 
to the mixed relative conditions downstream of the rotor. The output consists of rotor 
exit free-stream Mach number (XMFS1), rotor exit spacing (SPACING) in m (ft), 
trailing-edge thickness (TE) in m (ft), Mach number (XM2), axial Mach numbers (XMX1 
and XMX2), critical velocity ratios (V/VCR1 and V/VCR2), flow angles measured from 
axial direction (ALPH1 and ALPH2) in degrees, total -to -total pressure ratio (PT2/PT0), 
total-to-static pressure ratio (PT0/P2), static -to-total temperature ratio (T2/TT0), 
rotor kinetic energy loss coefficient (EBAR2), and rotor efficiency (ETA-N). Also 
printed are some miscellaneous parameters for the corrected blade profile including the 
final outlet flow angle and solidity. 


PROGRAM DESCRIPTION 

The program SSROTR designs supersonic rotor blades corrected for boundary -layer 
displacement thickness. The program uses the input value of the outlet flow angle 
(which controls the outlet circular turning) for the initial design. If the corrected outlet 
spacing is not equal to the inlet spacing, the program estimates a new value of outlet 
flow angle and repeats the design. The iteration procedure is essentially a half-interval 
search. The program input consists essentially of the surface Mach number distribution 


15 


and total flow conditions. The primary output is the corrected blade profile and the 
boundary-layer parameters. 

The program SSROTR is logically divided into two parts: (1) an ideal rotor design 
section, and (2) a section to calculate the boundary-layer parameters and aftermixing 
conditions. The ideal rotor design section consists principally of subroutines ROTORU 
and ROTORR. The boundary-layer section is controlled by subroutine BLAYR. The 
major subroutines called by BLAYR are PRECAL, LAMNAR, TURBLN, and AFMIX. 


SUBROUTINE ROTORU 


Subroutine ROTORU calculates the upper- and lower-surface transition arcs by the 
method of characteristics. The description of the method and the main dictionary of 
variables are given in reference 2. This routine is called only once for any given blade 
design. The program variables for ROTORU that are not in reference 2 are 


ETALOW 

ETAUP 

GAMMI 

GAMPI 

STML 

ST MU 


angle of tangent to lower surface, measured from axial direction 
angle of tangent to upper surface, measured from axial direction 
GAM-1, y - 1 
GAM+1, y + 1 

lower-surface critical velocity ratio 
upper-surface critical velocity ratio 


SUBROUTINE ROTORR 

Subroutine ROTORR calculates the rotated ideal passage coordinates, as well as the 
inlet and outlet passage spacing based on the current value of outlet flow angle. This 
routine is called once for each iteration of outlet flow angle. The description of the 
method and the main dictionary of variables are given in reference 2. The variables 
for ROTORR that are not in reference 2 are 

ETACL angle of flow along lower-surface circular arc 
ETACU angle of flow along upper-surface circular arc 
ETALN angle of flow along lower-surface inlet transition arc 
ETALO angle of flow along lower-surface outlet transition arc 
ETASN absolute value of inlet flow angle 
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ETASO 

ETAUN 

ETAUO 

SMLN 

SMLO 

SMUN 

SMUO 


absolute value of outlet flow angle 
angle of flow along upper-surface inlet transition arc 
angle of flow along upper-surface outlet transition arc 
lower-surface critical velocity ratio for inlet 
lower-surface critical velocity ratio for outlet 
upper- surface critical velocity ratio for inlet 
upper -surface critical velocity ratio for outlet 


SUBROUTINE BLAYR 

Subroutine BLAYR provides calls for the calculation of the boundary -layer charac- 
teristics and the mixing losses. The major subroutines called by BLAYR are PRECAL, 
LAMNAR, TURBLN, and AFMIX. These, in turn, call several other subroutines. The 
methods used for calculating the boundary-layer characteristics and the program de- 
scription are given in reference 4. The Lagrangian interpolation routine of reference 4 
has been replaced by a linear interpolation. 

Since a highly favorable pressure gradient exists for some portions of the blade 
surfaces, the range of the equations of reference 6 were extended by the method given 
in reference 6. These changes were made in LAMNAR and are for SW (temperature 
function at the wall) = 0. These changes are noted in the program with comment cards. 
Two curve-fit ranges were also extended. RCRIT and DIFF were extended as follows: 
RCRIT = 8. 3163 when SHAPE is greater than 0. 07; and DIFF = 44 000 KBAR + 700 
when KBAR is greater than 0. 03. 


SUBROUTINE AFMIX 

Subroutine AFMIX takes the boundary-layer parameters (displacement and momen- 
tum thicknesses) and the free-stream conditions at the rotor exit and calculates the 
aftermixing conditions by the method described in reference 5. In this loss model the 
flow sufficiently downstream of the blades is assumed to be mixed to uniform conditions. 

The calling sequence for AFMIX is CALL AFMIX (ALPHI, DELS, DELP, THETS, 
THETP, TE, SP, XMFSI), where 

ALPHI rotor exit flow angle measured from axial direction 

DELS displacement thickness at rotor exit for suction (upper) surface 

DELP displacement thickness at rotor exit for pressure (lower) surface 
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THETS momentum thickness at rotor exit for suction surface 

THETP momentum thickness at rotor exit for pressure surface 
TE trailing-edge thickness 

SP exit spacing 

XMFSI free-stream Mach number at rotor exit 

The subroutine description and the main dictionary of variables are given in reference 3. 


PROGRAM LISTING 


$ I EFTC SSROTR DECK 

C COMBINED BOTCP- BOUND ARY LAYER P PC CP AM 
CO MMCN/COD E/KCP E, DELS RL, DELS RU 

CO.MMON/CTOBL/ALPH1 ,SPA , XMAX , m E r STUB BU , NT U BPI. , CT HET U , CT H ETL 
COMMON/ LNK1/XOML (100) ,YOMl (100) ,S M L(100) ,ETAML(100) ,NI,,XOMU (100) , 
1 YOMIJ ( 1 0 0) , SMU (1 00) , ETA MU (100) , NO 

COMMON/LNK2/SMIN, SMO UT , R ET A T, CST A P , S TG M A O , G ST APO , G ST A PI , G A M M A , 

* CONVER , FECCNV 

COMMON/ FGT1/V (4) , E ET A N , N K J ( 0 ) , D R (4) 

COMMON/C 1/DU Ml (7) , MT UP B , DUM 2 ( 6) , CT HET , DU M 3 ( 0 1 0) 

DIMENSION X (100) , Y (100) 

1 BETAP = 0. 

BET AM ■= 0. 

C CALCULATE UNFOTATED ELADE COORDTNATFS 
CALI. POTOPU 
CALL INPUT1 
CALB = 0. 

C CALCULATE ROTATED BL.ADF COORDINATES 

2 CALL PCTCRR (CALB) 

CALB = 1. 

C NORMALIZE ROTATED PLAPE COORDINATES TO CHORD LENGTH 
DO 0 1=1, NO 
X(I) = XC.MU (I) /CSTAF. 

3 Y (I) = YCMU (I) /CS TAP 

KODE = 0 

NTURB = NTUPPU 
CTHFT = CTHETU 

CALL INPUT? ( X , Y , S M U , ET A MU , NU ) 

C CALCULATE UPPER SURFACE EOUND^FY LAYF? 

CALL PLAYR 

C NORMALIZE ROTATED BLADE COORPINATFS TO CHORD LENGTH 
DO 4 1=1, NL 
X{I) = XO.ML (I) /CSTAR 
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4 Y ft) = YCML (I) /CSTAP 
K01E = 1 

NTTIPB = NTURBL 
CTHET = CTHETL 

CALL TNPIIT2 (X,Y,SML, ETAML,NL) 

C CALCULATE LOWEF SURFACE BOUNDARY LAYER 
CALL ELAYR 

C CHECK FOR EQUAL SPACING AT INLET AND OUTLFT 
DELTOT = (DELSRL + DELSPU) *CSTAR 
OUTSPA = GSTAPO + DELTOT/CO S ( E ET AT*CON VEP) 

ELIMIT - . 0C01*CSTAR/XM AX 

IF (ABS (CST ARI - OUTSPA) .LF. ELIMIT) GO TO 5 
IF (OUTSPA .GE. GSTARI) EETAP = E E' p AT 
IF (OUTSPA .LT. GSTARI) BETAM = BET AT 
GSTAPO = OUTSPA 

IF (EETAP „NE. 0. .AND. BETAM .ME. 0.) GO TO 6 
C CALCULATE NFW OUTLET FLOW ANGLE 

BETAT = — ARCOS (COS (BETAN) *SMIN/SMOUT* ( ( 1 . + (GAMMA— 1. ) /2. * 

* SMOUT*SMOUT) /{I. + (GAMMA-1.) /2.*SMIN*SMTN) ) ** ( (GAMMA+ 1.) / (2.* ( 

* GAMMA-1.))) + DELTOT/GSTAPT) *PFCONV 

GO TO 2 

6 BETAT = (BETAP + BETAM) /2. 

GO TO 2 

5 WRITE (6,100) BETAT, GST API , OUTSPA 

100 FORMAT (//5X, 19HITER ATION COM PL ETED/7 X , 1 1 H 3 ET A (OUT) =,F10.4,4H DEG 
*, 10X, 8HG* (IN) =,F 1 0. 5, 1 C X , 22HG* (OUT) + DEL (TOT AL) =,F10.5) 

CALB = 2. 

CALL POTORP (CALB) 

GO TO 1 

END 

$ I E FTC INPU1 DECK 

SUBROUTINE INPUT 1 

C POUNDAPY LAYEP INPUT 

COMMON/LTU/LD AT 

COMMON/L NK2/UPKAC, SMOOT, A LP1M, CSTAP, ST G MAO, GSTAPO, GST a PI, GAMMA, 

* CCNVFR, FECGNV 

COMMON/CTOBL/ALPH1 ,SPA, X M A X , T E , NT U P BU , NT U FB L , CTH ETU , CTH ETL 
COM MCN/C 1/GAM, R , PTZ, TT7. , (JPK ACH , N ST , N VP , NTURB , KPVM, KEM , KSMTH , 
1KSPLN,KLF,KATCH,CTHFT, DLAM,TLAK , DT UP p , TTU PB , KPP E , KGP AD , KS DF , KL AM , 
2KMAIN,KPROF,X (100) , Y (100) , PPFS ( 100) , UF (1 00) , ME (100) , POPTZ (100) , 
3VOVCP (ICO) , TWAL (100) , ETA (10O) 

READ (LDAT,1P0) NTUPBU, NTU P B L , M V P , F , PT 7 , TTZ ,XMAX,TE, CTHE^U , CTHETL 
100 FORMAT (313, IX, 7 (F9.5, IX) ) 

READ (LDAT , 1010) KPP E , KGP AD , KS D F , K I. AM , K M ATM , K PROF, K EM 
1010 FORMAT (715) 

GAM = GAMMA 
UPKACH = UPMAC 

RETURN 

END 
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SIEFTC ROTO LIST, DECK 
SUBROUTINE EOTORU 

C**** CALCULATE THE UNROTATED BLADE COORDINATES 
COMHON/LTU/LDAT 

CO NEON/ FACT OR/PERM , GAMM 1,GAMP1,GAM 

CONNON/ROT1/VIN, VOUT, VLCW , VUP , B ET A N , NPEP, KKAXN, KM A XO, JN AXN , t7M AXO , 

1 PIN, ROUT, FLOW, RUP 

CON MON/ROT 2/XLOW (400) ,YLOW (40 0) , ST EL (400) , FT ALOW (400) , K NDEX , 

1 XU P (40 0) , Y UP (4 00) , STM 0 (400) ,ETAUP (400) , J N D E X 
COMMON/LNK2/SMIN,SMOUT,PETAT,CSTAR,SIGMAO,GSTARO,GSTARI, GANNA, 

* CONVER,RECCNV 

LOGICAL ANGLE 

EXTERNAL FOFRS 

DATA LD AT/5/ 

F (V ,FN) - 2. * V - HALFPI* (PERM- 1.) - 2. * (FN- 1 . ) *DELV 

1 READ (LDAT ,11) BET AN , VI N , VLOW , VUP , VOUT , B ET AT , DELV , G AM 
11 FORMAT (8 (F6. 2, 2X) ) 

CC CONVERSION FACTORS AND CONSTANTS 

GANNA = GAM 

CONVER = . 174532925F-0 1 
RECONV •= 57.2957796 
HALFPI = 3.14159265/2. 

C ONE POINT WILL BE PRINTED FOR EVERY NPER POINTS CALCULATED 

NPER = 10 

IF (DELV .GE. 0.2) NPER - 1 
GANP1 = (GAM + 1.)/2. 

GAMM1 = (GAM-1. )/2. 

PERM = SQRT (GAMP1/GAMM1) 

X0 = 1 . /PERM 

X 2 = 0. 999999999 

XT NTI. = (X0 + X2)/2. 

ANGLE = .TRUE. 

IF (VLOW .LE. AMIN1 (VUP, VIN, VOUT) ) GO TO 120 
WRITE (6,119) 

119 FORMAT (//3 IX, 70HV (LOW) MUST BE LESS than OR EQUAL TO THE MINIMUM 
IGF V (UP) ,V (IN) , V (OUT) ) 

ANGLE = .FALSE. 

120 IF (VUP .GE. AMAX1 (VIN , VOUT) ) GO TO 118 
WRITE (6,117) 

117 FORMAT (//33X,66HV (UP) MUST PF GREATER THAN OR EQUAL TO THE MAXIMU 
1 M OF V (IN) , V (OUT) ) 

ANGLE = .FALSE. 

118 VUMAX = HALFPI* (PEPM-1.) +RECONV 
IF (VUP .LE. VUMAX) GC TO 116 
WRITE (6,115) VUMAX 

115 FORMAT (//4 1 X, 37HV (UP) MUST BE LESS THAN V (UP) (MAX) = ,F9.4,4H DEG 
1) 

ANGLE = .FALSE. 
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116 IF (.NOT. ANGLE) GO TO 1 


WRITE (6,07) 

97 FORMAT (1H1,57X, 17 H DESIGN PARAMETFFS) 

CC MISCELLANEOUS CALCULATIONS 
DELV = DELV*CCNVEB 
FN = 1. 

V = V I N *CCN VER 
DO 4 1=1,2 
FOFX = F ( V , F N ) 

CALL FOOT (XO,X2,XINTL,FOFX,FOFRS,X1) 

IF (I .EO. 2) GO TO 4 
RIN = XI 

V = VOUT*CON VEB 
4 CONTINUE 

BOUT = XI 

SS MIN = 1./BIN 
SSMOUT = 1 . /BOUT 

SMS = SSMIN 
I = 1 

16 SM = SOFT ( ( (1./GAMP1) *SKS*SMS)/ (1.- ( G AMM 1/G AMP 1 ) *SMS*SMS) ) 

GO TO (17, 18, 19,20) , I 

17 SMIN = SM 
SMS = SSMOUT 
1 = 2 

GO TO 16 

18 SMOUT = SM 

DELV = DFLV*RFCONV 

CC TPINT ALL DESIGN PARAMETERS 

WRITE (6,95) BETAN,VIN,VUP,VOUT,UFTAT 
95 FORMAT (/2X, 10HBETA (IN) =,F8.4,4H D FG , 4 X , 7 H V ( IN ) =,FR.4,4H DEG,6X 

1 7 H V (UP) =, F9.4, 4 H DEG, 7X, 8HV (OUT) =,F8.4,4H DEG , 4 X , 1 1 HB FT A (OUT) 
2,F9.4,4H DEG) 

WRITE (6,94) DELV, VLOW, GAM 

94 FORMAT (/20X,9HDELTA V =,F8.4,4[i D EG , 1 1 X , 8 H V (LOW) = ,F8.4,4H DEG, 

1 1 1X,7UGAMMA = , F8. 4 ) 

CC CONVERT FROM DEGREES TO FABIANS 

VIN = VIN*CCNVFR 
VOUT = VOUT+CCNVEP 
VUP = V U P*C C N VER 
VLOW = VLOW *CON V F B 
BETAN = BETA N*CCNVER 
DFLV = D FLV *CCNV FF 

CC CHOOSE LONGFST TRANSITION ARC OF LOWER SURFACE 
VNL = VIN - VLOW 
KMAXK = (VNI/DEL V) + 0.6 
VOL = VOUT - VLOW 
FMAXO = (VOL/DELV) + 0.5 
KMN = MAXO (KMAXN, KMAXO) 

V - A M A X 1 (VTN, VOUT) 



CC CALCULATE q * (LOW) = R L 0 W , M * ( LC V) = SSM LC W , K (I 0 W ) =S Y. LO V 
IF (VLOW . EQ . 0.0) GO TO 2 
FN = KMN + 1 
FO^X = F ( V , F N ) 

CALL FOOT (XO,X2,XINTL,FO r X,FOFP.s f PLOW) 

GO TO 3 

2 BLOW = 1.0 

3 SSMLOW = 1./RLOW 
SMS = SSMLOW 
1=3 

GO TO 16 
19 SMLOW = SM 

CC SET INITIAL FOINTS FCR LCWFP APC C A I C UL A T ION S 
KNDEX = KMN/NPEP 
KDEX = KNDEX 
STML (KDEX+1 ) = SSMLOW 
XLOW (KDEX+1) = 0.0 
YLOW ( KDEX + 1 ) = PLOW 

ETALOW (KDEX+ 1) = 0. 

PHIKP1 = - (V-VLOW) + FLOAT (KMN) *DELV 

UMKP1 = AFSIN (SQRT (G AMP 1* RLOW*RLOW - GAMMl)) 

TXLO = XLOW (KDEX+1) 

TYLO = YLOW (KDEX+1) 

CC CHOOSE LONGEST TRANSITION ARC OF UPPER SURFACE 
VUT = VUP - VC UT 
JMAXO = (VUT/DET.V) +0.5 
VUI = VUP - VIN 
JMAXN = (VIII/DELV) +0.5 
JBN = MAXO (JMAXO, JMAXN) 

V ■= A MINI ( VCUT , V IN) 

CC CALCULATE P*(UP)=PnP, M * ( UP) =S S M UP , K(UP)=.MTJP 
FN = - (JMN+ 1) + 2 
FOFX = F (V, FN) 

CALL FOOT (XO , X2 , XINTL, FOFX ,FOFFS, PUP) 

SSHUP = 1./BUP 
SMS = SSMUP 
i = a 
GO TO 16 
20 SMUP = SM 

CC SET INITIAL POINTS FOR UPPER APC CAICIJLATIONS 
JNDEX = JMN/NPEP 
JDEX = JNDEX 
STMU (JDEX+1) = SSMUP 
XUP (JDEX + 1) = 0.0 
YUP ( JDEX+ 1 ) = PUP 

ETAUP (JDEX+1) = 0. 

PHI JP 1 - -(VUP-V) + FLOAT (JMN) *DELV 
UMJP1 = AFSIN (SORT (GAHP1*PUP*PUP - GAMMl) ) 

TXUP = XUP (JDEX+1) 

TYUP = YUP (JDEX+1) 

IF (VIN . EG- VLOW .AND. VLOW .EG. VOUT) GO TO HO 
C****CALCULATE COORDINATES FOP LOWER TRANSITION APC - UNROTATED 
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KDEX = KNDEX + 1 

NUM = 0 

V = A MAXI (VI N, VOUT) 

KEEP = KMN + 1 

DO 10 KK- 1 , K M N 

K = KEEP - KK 

NtJK = NUM + 1 

PHIK = PHI K E 1 - DEL V 

FN = K 

FOFX = F(V,FN) 

CALL ROOT (XO, X2 r XINTL, FOFX , FOFFS , TR) 

TX = TR*SIN (PHIK) 

TY = TR*COS (PH IK) 

EM WK = TAN (-PHIKP1) 

UMK = ARSIN (SORT (GAMP1*TR*TR - GAMM 1) ) 

EM K = -TAN ( (PHIK + UMK + PHIKP1+UMKP1)/2.) 

TEMP = TYLO - EMWK*TXLO 

TEMPP = TY - EMK*TX 

TEMPPP = EMK - EMWK 

TXLO = (TEMP - TEMPP) /TEMPPP 

TYLC = ( (EMK*TSMP) - f FM W K*T EMP P) ) /TEMPPP 

PHIKP1 = PHIK 

UMKP1 = UMK 

CC SAVE EVERY "NPER-TH" POINT 

N = NUM - (NOM/NPER) *NPEP 
IF (N .GT. 0) GO TO 30 
KDEX = KDEX - 1 
STML (KDEX) = 1./TR 
XLON(KDEX) = TXLO 
YLOW(KDEX) = TYLO 
ETALOW (KDEX) = -PHIK 

30 CONTINUE 

100 IF (VIN . EQ. VtIP .AND. VHP .EO. VOUT) GO TO 200 

C****CALCULATE COORDINATES FOP UPPER TRANSITION APC - UNROTATED 
JDEX = JNDEX + 1 
NUM = 0 

V = A MINI (VOUT, VIN) 

JFEP = JMN + 1 

DO 41 JJ=1,JMN 

J = JEEP - JJ 

NUM = NUM + 1 

PHIJ = PHIJF1 - DEL V 

FN = -J + 2 

FOFX = F (V, FN) 

CALI ROOT ( X 0 , X2 , XINTL, FOFX ,FOFRS , T R) 

TX = TR*S I N (FHIJ) 

TY = TR*COS (PHIJ) 

EMWJ = TAN (-PH1JP1) 

UMJ = APSTN (SOFT (GAMP1*T n *TF - G A M M 1 ) ) 

EM J = TAN ( (— PHIJ + UMJ— PHI JP1 + UMJP1) /2.) 

TF"P = TYUP - E M W J * T X II P 

TEMPP = TY - E M J * T X 

TEMPPP = EMJ - FMWJ 

TXUF = (TEMP - TEMPP) /TEMPPP 

T Y U ^ = ( ( EM J*TEM P) - (EMU J*TFMPP) ) /TFMFPP 



PHI JP 1 = PHI J 
UMJP1 = UK J 

CC SAVE EVFFY "NPER-TH" POINT 

N = NIJK - (NUK/NPER) *NPEF 
IF (N .GT. 0) GO TO 41 
JDEX = JDEX - 1 
STMU(JDEX) = 1./TR 
XUP(JDEX) = TXUP 
YUP (JDEX) = TYUP 
ETA UP (JDEX) = - PHI J 

41 CONTINUE 

CC MISCELLANEOUS OUTPUT 
200 WRITE (6,622) 

622 FORMAT (//54X, 24HE ISCELL A N EOUS PARAMETERS//) 

WRITE (6,1000) SSMTN,SMTN, SMOOT, SSMCUT 

1000 FORMAT (/25X,8HM* (IN) = , F9 . 4 , 3 X ,7 HK ( IN) = ,F9. 4, 10X, 8HM (OUT) =,F9.4 
1,5X,9HM* (OUT) =,F9.4) 

WRITE (6 , 1 CC 1 ) SSMLCW,S MLOW ,SMUP,SSMUF 

1001 FORMAT (/25X, 9HM* (LOW) =, F9 . 4 , 2X , 0 HM (LOW) =, F9. 4 , 1 1 X, 7HM (UP) =, 

* F9.4,6X,8HK* (UP) =,F9.4) 

KNDEX = KNDEX + 1 
JNDEX = JNDEX + 1 

RETURN 
END 

$ I EFTC ROO DECK 

SUBROUTINE BOOT (X0 , X2, XINTL, FOFX ,FUNC, X 1) 

DOUBLE PRECISION X,XX0,XX2 

C WE ARF SEEKING AN X SUCH THAT FUNC(X) = FOFX WHERE 
C FUNCTIONAL VALUE 

C 1 LOCATE FOFX IN (F0,FX) OR (FX,F2) WHERE FX IS 

C APPROXIMATION TO FOFX 

C 2 LET X = 1/2 (XX0 + X) OR X = 1/2(X + XX2) 

C 3 IS FUNC(X) = FOFX ? IF NOT, REPEAT PROCEDURE 

XXO = XO 
XX 2 = X 2 
FO = EUNC (XXO) 

F2 = FUNC (XX2) 

IF ( FOFX .LT. FO .AND. FOFX .IT. F2 .OP. FOFX .GT. FO .AND. 

1 FOFX -GT. F 2 ) GO TO 1005 

IF (ABS (FOFX-FO) .LF. .0000 1) GO TO 1007 
IF (ABS (FOFX-F2) .LE. .00001) GO TO 1008 

X = XINTL 
KOUNT = 0 
1C00 XI = X 

KOUNT = KOUNT + 1 
A = FOFX - F2 
FX = FUNC (X) 

IF (KOUNT .GE. 60) KFITF (6,1004) KOUNT , X , FX , FOFX 


FOFX IS A KNOWN 
THE PREVIOUS 


24 



n n 


1 HO 4 FO R.M AT ( 1 H L , 9 H KOUNT ,G16.9,9H X ,0 16.0,911 FX ,G16.9, 
1 9 H FO F X , G 1 6 . 9 ) 

IF (AES (FX-FOFX) .LE. .00901) RETURN 
IF (KCONT .E0- 75) GO TO 1002 
IF (A* (FX-FCFX) .LT. 0.) GO TO 1001 
X X 0 = X 

X = (X+XX2) /2. 

GO TO 1000 
1C01 XX2 = X 

X = (XXO+X) /2. 

F2 = FX 
GO TO 1000 

1002 WRITE (6,1003) 

1003 FORMAT (//30X,62H75 ITERATIONS HAVE BEEN PERFORMED WITHOUT CONVERG 
1ING TO A ROOT) 

RETURN 

1 C05 WRITE (6,1006) FOFX 

1006 FORMAT (// 1 OX , 7 H F ( X) = ,G16.9,31H IS OUTSIDE OF SPECIFIED LIMITS) 

RETURN 

1 C07 XI = XO 
RETURN 

1008 XI - X 2 
RETURN 

END 

$ I BFTC FELI DECK 

FUNCTION FOFRS (X) 

DOUBLE PRECISION X 

COM MON/FACTOR/PERM, GAM Ml ,GA KP1 , GAM 

AEG1 = 2.*GAMM1/(X*X) - GAM 
ARG2 = 2. *G A MP 1 * X * X - GAM 

IF (ABS(ARGI) .GT. 1.0 .OR. AES ( A RG2) . GT . 1.0) WRITE (6,1) APG1 

1 , A R G 2 

1 FORMAT (// 14X,61HARGUMEKT OF ARCS IN IS OUTSIDE DOMAIN OF DEFINITIO 
IN ARG1 = , G 1 6 . 9 , 11H ARG2 = ,G16.°) 

FOFRS = PERM*ARSIN (ARG1) + AFSIN(AFG2) 

RETURN 

END 

$ I PFTC ROTP LIST, DECK 

SUBROUTINE DOTORR (CALB) 

CALCULATE (1) ROTATED BLADE COORDINATES 

(2) INLET AND OUTLET ELATE SPACING 

DIMENSION XLOWN (80) , YLOWN (80) ,SMLN (80) , ETALN (80) ,XCI,OW (35) , 

1 YCLOW (35) , ETACL (35) ,XLOWO (80) , YLOKO (8 0) , SMLO (80) , ETALO (80) 
DIMENSION XSN (1 1) , YSN (1 1) , ETASN ( 11) , XUFN (80) , YU PM (80) , SHUN (80) , 

1 ET At! N (80) , XCU P ( 35) , YCUP (35) , ET ACU(3 5) ,XUPC (80) ,YUPO(R0) ,SMUO(8C) , 

2 ETAUO (8 0) , XSO (1 1) , YSC ( 1 1) , FTASO ( 1 1) 

COMMON/LNK1/XOMT. ( 1 00 ) , Y 0 ML ( 1 00) , S ML ( 1 0 0) , FT A ML ( 1 00 ) , NL , XOMU ( 1 00) , 


1 YOMU (ion) ,SKU ( 1 0 O') , ETA PM (ICO) , MU 
COMKOK/POT1/VTK, VOUT r VLnK,VtJP,nFTAP r KPEP , KM A XN, V M A XO, Jfl A XN, JM A XO, 

1 KIN, FOUT, PLOW, RUP 

COMMCN/P nTP/XLOlf («D0) , YI.OW (4^0) ,STM I. (4^0) , FT A 10 W (400) , K N D E X , 

1 XUP (400) , YtJP (400) , S T M U (40n) , FT A HP (40P) ,.tsDEX 

C0MM0N/LKK2/SMTN,SM0UT, P F T AT, CST" R , S TGM AO , 0 ST A RO , G ST A p I , 0 A M « A , 

* CONVFR,FECCNV 

IF (CAL? .ME. 0.) GO TO 2 
WPITE (6,100) 3FTAT 

100 FORMAT (//10X, 15HINPUT BFTA(O) =,"10. 4, an Di'G///) 

GO TO 3 

2 IF (CALB .NE. 1.) GO TO 3 
WRITE (6,101) SETAT,r,STAPO 

101 FORMAT (// 10X,20HCALCUI. ATED B FT A ( 0) =,F10.4,4K DFG , 5X , 1 4HWH EN G* (0 
HIT) = , F 1 0. 6) 

3 BETAT = BETAT*CONVER 


ALPHLN = (YIN - VLOW) - PETAL’ 

ALPHLO = -(VOUT - VLOW) - B FT A T 

IF (ALPHLN .LE. 0. .AMD. ALPHLO .GE. 0.) GO TO 4 
WRITE (6,102) 

102 FORMAT (//27X , 79HV (LOW) HOST PE GRFATFR THAN OP FPIJAL TO V(IN) - I? 
1 E T A (IN) AMD V (OUT) + BETA (OUT)) 

STOP 


C****CALCUI.ATE 
4 SINAI , V - 
COS ALN = 
ABSALN = 
SINALO = 
COSALO = 
ARSALO = 


COCPDI NATES 
SIN (ALPHLN) 
COS (ALPHLN) 
ABS (ALPHLN) 
SIN (ALPHLO) 
COS (ALPHLO) 
ABS (ALPHLO) 


KDEX = KNDEX + 1 


FCF LOWFP TRANSITION 


KN = (KKAXN/NPFR) + 2 
KO = (KM A XO/NPFP) + 2 


A PC 


ROTATED 


DO 6 KK=1,KNDFX 
K = KDEX - KK 
KN = KN - 1 
KO = KO - 1 

IF (KM .LE. 0) GO TO 6 
SMLN(KN) = STML(K) 

XLOKN(KN) = YI.OW (K) *STNAI N + XI CW ( K ) * COS » IN 

YLOWN(KN) = YLOW (K) *rOSAL" - X I CW ( K) * S I N A L N 

ETA.LN(KN) = FT A LOW (K) + ABS AIN 

6 IF (KO .IE. 0) GO TO 5 
SHLO(KO) = ST ML (K) 

XLOWC(KO) = YLOW (K) * S I N A L 0 - X L CW ( K) *Cn$ A LO 

YLOWO(KO) = YLOW (K) *COS AID + XL CW ( K) *S I N AI.O 

FT A LO (KO) = ETALOW(K) + A55ALO 
6 CONTINUE 


ALPHUN = (VUP - VIN) - BET AN 
ALPHUC = -(VUP - VOP'P) - BETA.' 7 ’ 

IF (ALPHUN .LE. 0. .AMP. AIPHUO .GE. 0.) GO TO 7 
WRITE (6,10 3) 
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103 FORMAT (//2X, 75HV (UF) MUST RE LESS THAN OP KCUAL TO V(IN) + FETA (I 
IN) AND V (OUT) - RFTA (OUT) ) 

STOP 

C*+**CALCULATE COORDINATES FOP UPPER TRANSITION ARC - ROTATED 

7 SINAUN = SIN (ALPHUN) 

COSAUN = COS (ALPHUN) 

ABSAUN = ABS (ALPHUN) 

SINAUO = SIN (ALPHUO) 

COSAUO = COS (ALPHUO) 

ABSAUO = ABS (ALPHUO) 

JDEX = JNDEX + 1 

JN = (JMAXN/NPER) + 2 

JO = (JMAXO/NFEP) + 2 

DO 8 JJ=1, JNDEX 

J = JD EX - JJ 
JO = JO - 1 

JN = JN - 1 

IF (JO -LE. 0) GO TO 9 
S MUO (JO) = ST M U ( J) 

XUPO(JO) = YUP (J) * S I N A U C - X UP ( J) *COS A UO 
YUPO(JO) = YUP ( J) *CCS AUC + XUP ( J) *S I NAUO 
ETAUO(JO) = FTAUP(J) + ABSAUO 
9 IF (JN -LE. 0) GO TO 8 
S MU N ( J N ) = ST.MU (J1 

XUPN(JN) = YtJP (J) *SI NAUN + XUP ( J) *COS AUN 
YUPN(JN) = YUP ( J) *COSAtJN - XU P ( J) * S I N A UN 
ETAUN(JN) = ETA UP ( J) + ABSAUN 

8 CONTINUE 

C****CALCULATE THE INLET AND OUTLET ( D I M F NSIO N L ES S) BLADE SPACING 
YLASTI = YU PN f 1 ) ♦ TAN (EFT AN) * (XLOfc'N (1) - XHPN(D) 

GSTARI = YLCWN(I) - YLASTI 

YLASTO = YUPO ( 1 ) + T AN ( BET AT) * ( XLOKO ( 1 ) - XTTPO(I)) 

GST ARO = YLCWO(I) - YLASTO 

DALPH = 5. *CCNVER 
C ****CIRCULAR ARC (LOWER) 

ALPH = ALPHLO + DALPH 
ALPLOW = ALFHLN 
KOUNT = 0 

10 KOHNT = KOUNT + 1 

XCLOW (KOUNT) = FLOW*SIN (ALPLOW) 

YCLOW (KOUNT) = PLOW*COS (ALPLOW) 

ETACL (KOUNT) = ABS (ALPLOW) 

ALPLOW = ALPLCW DALPH 

IF (AES (ALPH- ALPLOW) . LF. .001) GO TO 11 

IF (ALPHLO -IT. ALPLOW -AND. ALPLOW . LT. ALPH) ALPLOW = ALPHLO 
GO TO 10 

11 NCL = KOUNT 

C****CIRCULAR ARC (UPPER) 

ALPH = A I PH UC * DALPH 
ALPHtJP = ALFHUN 
KOUNT = 0 
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12 KOUNT = KOUNT + 1 

XCUP (KOUNT) = RUP*STN (ALPHUP) 

YCUP (KOUNT) = RIJ P*COS (ALP HUP) 

ETACU (KOTJNT) = ABS (ALPHUP) 

ALPHUP = ALPHUP + DALPH 

IF (ABS (ALPH-ALPHUP) .LE. .001) GO TO 13 

IF (ALPHUO .LT. ALPHUP .AND. ALPHUP .LT. AIPH) ALPHUP = ALPHUO. 
GO TO 12 

13 NCU = KOUNT 

C****CALCULATE COORDINATES FOP STRAIGHT LINE PORTION OF UPPER ARC 
KOUNT = 1 

IF (X LOW N { 1 ) .LE. XUPN ( 1 ) ) GO TO 15 
WRITE (6,106) 

106 FORMAT (8H0 (FOTRP) , UX, 7 1 HUPPER SURFACE TNLFT LONGER THAN LOWER SUR 
♦FACE INLET - CASE TERMINATED) 

STOP 

15 DELXT = (XUPN ( 1 ) - XLOWN (1) ) / 1 1 . 

IE (XLOWO(I) .GE. XUPO(1)) GO TO 16 
WRITE (6,105) 

105 FORMAT (8H0 ( ROTRE) , 4 X , 7 3 HUPPER SURFACE OUTLET LONGER THAN LOWER SU 
*R FACE OUTLET - CASE TERMINATED) 

STOP 

16 DELXO = (XLCWO(I) - XUPO(1))/10. 

XSIN = XUPN (1) 

YSIN = YUPN (1) 

XSOUT = XUPC(1) 

YSOUT = YUPC ( 1 ) 

XSN (KOUNT) = XSIN 
YSN (KOUNT) = YSIN 
ETASN (KOUNT) = ETAUN(I) 

XSO (KOUNT) = XSOUT 
YSO (KOUNT) = YSOUT 
ETASO (KOUNT) = ETAUO(I) 

TANBN = TAN (EETAN) 

ABSBN = AES (EETAN) 

TANBO = TAN (BETAT) 

ABSPT = APS (BETAT) 

DO 14 1=2,11 

XSIN = XSIN - DELXI 

XSOUT = XSOUT + DELXO 

YSIN = YUPN ( 1 ) + TAM BN* ( XS IN - XUPNfl)} 

YSOUT = YUPO { 1 ) + TANPO* (XSCPT - XUPO(D) 

ETASN (I) = ABSBN 

XSN ( I ) = XSIN 

YSN (I) = YSIN 

ETASO (I) = ABSBT 

XSO(I) = XSOUT 

14 YSO (I ) = YSOUT 

C****FREPARF DATA FOR BOUNDARY LAYER PROGRAM 

CSTAR = SORT ( (XLOWO ( 1) - XLCWN(1))**2 + (YI.OKO(I) - YLO W m ( 1 ) ) ** 2 ) 

SIGN A I = CSTAR /GST APT 
SIGMAC = CSTAF/GSTAFO 

WFITE (6,104) GSTARI,STGMAI, CSTAR, GST AF 0 , S I GM AO 
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104 FORMAT (//39X, PHG* (IN) = , ?9 . 4 , 3 7X , 1 1 H SIR K A ( I N) = , F9 . 4 /6 9 X , 4 HC* = , 
* F9.4/39X,9HG* (OUT) = , F9. « ,36 X , 1 2HS IGMA (CUT) = ,F9.4//) 

13 FT AT = EET AT*RECONV 

SSMLOV = 1./RLOW 

SSMUP = 1./BCIP 

SSMIN = 1 . /PIN 

SSMOTJT = 1./RCUT 

IF (CALB -EC. 2.) RETOFN 

C STORE ROTATED ELADE COORDINATES FOR LOWFR SURFACE 
KN = (KMAXM/NPER) + 1 
DO 310 1= 1 , K N 
XOML (I) = XLOWN (I) 

YO ML (I) = Y LOW N (I) 

SML(I) = SHI.N (I) 

ETAML(I) = ETALN(I) 

310 CONTINUE 
ML = KN 
NCL - NCL-1 

no 311 1=2 f NCL 
NL = NL + 1 
XOML(NL) = XCLOW(I) 

YOML(NL) = YCLOW(T) 

SML(NL) = SSMLOW 
ET AML (NL) = ETACL(T) 

311 CO NT I NU F 

KO = (KMAXO/NPF.R) + 1 

DO 312 1=1, KO 
J = KO + 1 - I 
NL = NL + 1 

XOML(NL) = XLOWO(J) 

YOML(NL) = YLOWO(J) 

S M L (ML) = SMLO(J) 

ET AM L (NL) = ETALO(J) 

312 CONTINUE 

r STORE ROTATET ELADE COORDINATES FOP UFFFP SURFACE 
DO 313 1=1,10 
J = 12-1 
XOMU(I) = XSN(J) 

YOMU(T) = YSN(J) 

SMU(I) = SSMIN 
ET AMU (I) = ETASN(J) 

31 3 CONTINUE 
NU = 10 

JN = (JMAXN/NPER) *■ 1 
DO 314 1= 1 , J N 
NU = NU + 1 
XOMU(NU) = XUFN(I) 

YOMU(NU) = YUPN(I) 

S MU (NU) = SMTJN(I) 

ETAMU (NU) = ETAUN(I) 

314 CONTINUE 
NCU = NCU-1 

DO 3 1 S 1=2, NCU 
NU = NU + 1 
XOMU(NU) = XCUP(I) 

YO M U ( N U) = YCUP(T) 

SMU(NtJ) = SSMUP 
ETAMU (NU) = FT AC U (I ) 
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315 CO NT I MtJ F 

JO = (JMAXO/NPER) + 1 
DO 316 1=1, JO 
J = JC + 1 - I 
NU = NU + 1 
XOHU(NU) = XUPO(J) 

YOMU(NU) = YUPO(J) 

SMU(NU) = SHUO(J) 

ETAMU(NU) = FTAUO(J) 

316 CONTINUE 

DO 317 1=2 , 11 
NU = NU + 1 
XOMU(NU) = XSO(I) 

YOMU(NU) = YSO(I) 

SMU(NU) = SSMOUT 
ETAMU(NU) = ETASO(I) 

317 CONTINUE 

RETURN 

END 

STEFTC INPU2 DECK 

SUBROUTINE INPUT2 (X M , YM r V VCR , FT A H , NN) 

LOGICAL ERROR, TRANS, SEPFN 
CON NON/COD F/KCDE, DELS RL,D EL SRU 

CONMON/L NK2/U PH AC , SNOUT, ALP1M,CSTAR, SIGH AO, GSTARO, GST AR I , GAM HA, 

* CONVER, FECCNV 

CO MNON/CTOBL/ALPH 1 ,SPA,XMAX,TF, NTURBU , NTU PB L , CTHETU , CTHF.TI. 

CON NON/C 1/GAM, R , PTZ , TTZ r UPM A C H , N ST , M VP , NT UP 3 , KPV« , KEN, KS NTH, 
1KSPLN,KLE,K ATCH ,CTHET ,DLAM,TLAN,DT U R E ,TTU R B , KF RE , KGR A D, KS DE , KL AH , 
2KMAIN,KPROF, X (100) , Y (100) ,PPES (100) ,UF { 100) , HE (100) ,POPT7, (100) , 
3VOVCP (100) , TW A L { 1 00) ,ETA{100) 

COMMON/C3/XOM (1 00 ) , YOH (1 00) , S (100) , SOL ( 100) , AE (100) ,TSE(100) , 
1TAWL (100) ,TAWT (100) ,TBAR (100) ,PW (ICO) r SW (100) ,SUTHL (100) , 

2RHSW (100) ,RHSE (100) , HEADW (100) , MEADE (100) , NOW (100) , NU BAP (100) , 

3AA (100) , R3 (100) , FF (100) , DUDS (100) , D H O S (100) ,DMDL (100) 

COMMON/C 9/EPROR, TRANS, SFPRN 

DIMENSION XM (NN) , YN (NN) , VVCP (NN) , ETAN (NN) 

THE V/VCR DISTRIBUTION IS PECEIVFD FRON THE POTCP PROGRAM 
IT WILL BE USED TO CALCULATE THE PPESSU^F DISTRIBUTION 

ERROR = .FALSF. 

TRANS = .FALSE. 

SEPRN = .FALSE. 

C BOUNDARY LAYER SETUP 
KATCH = 0 
KS PL N = 1 
KLE = 1 
TLAH = 0. 

DLAM = 0. 

TTUPB = 0. 

DTURB = 0. 

NST = NN 
DO 2 1=1, NST 
XOM (I) = XN (I) 
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YOM (I) - YM (T) 

X fT) = XOM (I) * XM AX 
Y ( I) = YOM(I)*XMAX 
THAT. (I) = TTZ 
ETA (I) = ETAM (I) 

2 VOVCP (I) = V VCR (I) 
ALPH1 = -ALE IK 
SPA = XKAX/SIGMAO 


IF (NODE . NE. 0) GO TO 3 
WRITE (6,100) 


100 

FORMAT (1H1,53 X,23H*** 
GO TO 4 

IJPPEP 

SURFACE 

***) 

3 

WRITE (6,101) 




10 1 

FORMAT (1H1,53X,23H*** 

LOWER 

SURFACE 



4 WRITE (6,1020) NST , N V P , C, A M , B , PTZ , TTZ ,TJ EH ACH , XK AX , NTURB, CT HET , TE 
1020 FORMAT ( 1 HO/SX , 22 HBOU N DA R Y LAYER - I NPUT///SX ,3HNST,5X, 3HNVP,7X, 

* 3HGAM, 1 1X , 1HR, 1 2X , 3H PTZ , 1 1 X , 3HTTZ , 0 X , 6H UPM ACH , 3 X, 4 HXM A X , 7 X , 

* 5HNTURB,6X , 6HCTH FTA ,7X ,2HTE/4 X ,2 (T3,5X) , F7 . 3 , 3 (5X , F9 - 2) , 5X,F8.4, 

* SX,FB.5,6X,I3,7X,F7.3,5X r F7.6///) 

WRITE (6,1090) KPRE,KGRAD,KSDF,KLAM, KKATN, KPROF , KEM 
1C9 0 FORMAT (/6X , 4HK PRE, 7X, 5 HKGR A D , 7X ,4 HKS DE , B X , 4HKL A M , 7X , 6H KM AI N , 7X, 6 HK 
1 PROF, 7X, 5H KEM /7X,l2,9X,I2,10X,r2,10X,I2,9X,l2,10X,l2, 10X,I2///) 
WRITE (6,1032) (XOM ( I) , YOM (I) , V VCR ( I) ,T W A L ( I) , 1= 1 , NST) 

1032 FORMAT (8X, 3KXOM,7X, 3H YOM , 9 X, 4 H V VCR ,9X,4HTWAL/(SX,F8.5,2X,F8.5,2X, 
1F1 1.5,3X,F9.2) ) 

IF (NST. GT. 1 CC.OR.NTURB.GT. NST. OF. KFM.IT. 0 . OP . K EM . GT. 1 - O R . K S PLN. LT . 
10. OR. KSPLN. GT. I.OR. KLE. LT.O.OP.KLF.GT.I.OP. K ATCH . LT . 0 . O F. K ATCH. GT. 
21) GO TO 70 
RETURN 

70 ERROR = .TRUE. 

WRITE (6, 1 17C) 

1170 FORMAT (////10X, 48HERROP IN INPUT DATA. FECHECK INPUT INSTRUCTIONS 

1) 

RETURN 

END 

SIEFTC BLAYFR DECK 

SUBROUTINE ELAYR 

C EOUNDARY LAYER MAIN PROGPAM 

LOGICAL ERRCP, TRANS, SERF N’ 

REAL ME 


COMMON/CODF/KCDE, DELSRL , DELSRU 

CO M MON/CTOBI./ AL , SP, XMAX,TE, NTU PB U , NTUR BL , CH TETU , CTH ETL 

COHMON/C1/GAM,P,PTZ,TTZ,UPM ACH, NST, NV F, NTU R E , KPV M , K EM , KS MT H , 
1KSPLN, KLF, NATCH, CTH FT, DL AM , TL A M , DT UR B , TTUPB , KPPE , KG R A D , KSDE , XL AM , 
2KMAIN,KPROF, X (1OO),Y(1OC),PRFS(1OO),HE(1OO),ME(1OO),POPTZ(1O0), 
3VOVCR (100) ,TWAL (100) , ETA (100) 

CO MMON/C 2 /PS Z , TSZ , UZ , AS 7 , AT 7 , R H SZ , REITZ , M USZ , M UTZ , NO SZ , N UT Z , CP , 
1PR,TC,ARCL 

COMMON/C3/XOM (100) ,YOH (100) ,S (100) , SOL (100) , AE(100) ,TSE (100) , 

IT AWL ( 100) ,TAWT (100) ,TBAP(100) ,RW ( 100) ,SW (100) ,SUTBL (100) , 

2RHSW (100) , PHSF (100) , HEADW (ICO) ,HFADE(100) ,NUK (100) , MUBAR ( 100) , 

3 A A (10 0) , BB (100) , FF (100) ,DUDS (100),DMDS(100) , DMDL MOP) 

COMMON/C4/THET (100) ,DELSR(100) , DELTA (ICO) , FORM (10 0) , 



1F0RMI (100) , FCRMTR ( 100) , PTH ( 100) ,R:1T 1(100), CF(100), 

2TAUK (100) , MUSS (100) , DTDY (100) ,KTRAN (100) ,CPN (100) 
CCMMCN/C9/FRROR, TRANS, SFPRN 

CALL PRECAL 
IF (ERRCP) RETURN 
CALL LAMNAP 
IF (ERROR) RETURN 
IF (SEPRN) GO TO 20 
IF (.NOT. TRANS) GO TO 20 
CALL TUPBLN 
IF (ERRCP) RETURN 
20 CALL PROFIL 

IF {KODE. EC. 1} GO TO 9 
DELS = DELSR (NST) 

THETS = THET (NST) 

KODE = 1 
RETURN 

9 DELF = DELSR (NST) 

THETP = THET (NST) 

ALPH1 = AL 
SPA = SP 

CALL AFMIX (A I, PH 1, DELS, DFLP, THETS, THETP, TE, SPA, ME (NST)) 

KODE = 0 
RETURN 

END 

$ T EFTC PRFCI LIST, DECK 
SUBROUTINE PRKCAL 

COMMON/ GAMP M/GMP1,GMM1 

COM MO N/CTOBL/DtJM (2) , XM AX , DUM2 (5) 

CO MKON/C 1/G A M , P , PTZ, TTZ, U PM ACH , NST, NVF, NTHRB, KPVM, KEM, KSMTH, 
1KSPLN,KLE,KATCH, CTHET, DI. AM , TL AM , DT UR B ,TTUPB , KPFE , KG R A D , KSDE,KLAM, 
2KMAIN,KPPOF,X (100) ,Y (100) ,PRES (100) ,UE(1 n 0) , ME (100) , PORT Z (100) , 
3VOVCR (100), TWAL (100), ETA (100) 

CO M MON/C 2/P SZ,TSZ,UZ, AS Z , ATT , RHSZ , PHTZ , M USZ , M UTZ , NUSZ , NUT Z , CP , 
1PR,TC, ARCL 

COMMON/C3/XOM (100) , YOM (100) , S (100) , SOL (100) , AE (100) ,TSE (100) , 

IT AWL (100) , TAWT (100) , TB A P (100) , PW (100) ,SW (IPO) ,SUTHL (10Q) , 

2RHSW (100) , RUSE (100) ,HEACW(100) ,HEADE(100) , NUW (100) , MUBAP ( 10 0) , 

3AA (100) , BB (100) , FF ( 100) , DUDS ( 100) ,DMDS(100) ,DMDL(100) 

COM MON/C 9/ERROR, TRANS, SEPRN 

DIMENSION SDER(IOO) , CMU (20) ,CPR(20) ,CTC (20) 

REAL KUSZ,MUTZ, NUSZ, MUTZ, MUSLE, MUS LM, HE, NUK , MTJFAP 
LOGICAL ERROR, TRANS, SEPPN 
C 

C READ DATA FOR MU, PR, AND TC CURVE FITS 
C 

DATA (CMU (I) ,1=1,5) /-.0194S170, 1.3 0196 31,-. 34511 323, 

1. 068277826,-. 00566593/ 

DATA (CP R(I), 1= 1,5)/. 8557,-. 234136, .1078624, 

1-. 0236214,. 00202863/ 

DATA (CTC( I) , 1=1, 5)/-. 0 38 39 3 2 3, 1.26 974 27, -. 3 00 112^2, 

1.08743781, -.00967472 6/ 

C 

GMP 1 = GAM + 1. 

GMM1 = GAM - 1. 
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INITIALIZE STATIC AND TOTAL PARAMETERS 

TSLE= S18.688 
TSLM= 2RR.160 
MUSLE- 3. 71 14 02E-7 
MOST. K= 1. 777029 E— 5 
TCSLE= 3. 202206E-3 
TCSL M= 2.561796E-2 
TSZ = TTZ/ (1 .+GMM1/2.*UPMACH**2) 

PSZ = PTZ* (TSZ/TTZ) ** (GAM/GMM1) 

RHSZ= PSZ/R/TSZ 
RHTZ= PTZ/P/TTZ 
ASZ = SQRT (GAM*H*TSZ) 

ATZ = SQRT (GAM*R*TTZ) 

UZ= U PM ACH* ASZ 
CP = F*GAM/GMM1 
IF (KEM. FQ. 1) GO TO 10 
TCON= 198.60 
TR1= TSZ/TSLE 
TP2= TTZ/TSLE 
GO TO 20 
10 TCON= 110.33 
TP 1= TSZ/TSLM 
T R 2= TTZ/TSLM 

20 CALL CUPVFT (CFP, PR, TR1 ,0,4,0) 

CALL CURVFT (CTC,TC,TR1, 0,4,0) 

CALL CUPVFT (CMO , M OSZ ,TR 1 ,0,4,0) 

CALL CUPVFT (CM U, M UTS, TR2, 0,4,0) 

IF (KEM . EQ. 1) . GO TO 30 
TC= TC*TCSLE 
MUSZ= MUSZ*MUSLE 
MUTZ= MUTZ* MUSLE 
GO TO 40 
30 TC= TC*TCSI.M 

MUSZ= MUSZ*MUSLM 
MUTZ= MUTZ* MUSLM 
40 N 1JSZ = MUSZ/PHSZ 
NUTZ= MUTZ/ P HT Z 

C CALCULATF GEOMETPY PATIOS AND A^C LENGTHS 
C 

S (1) = 0. 

DO SO 1=2, NST 

SO S(I)= S (1-1) +SQPT ( (X (I) — X (1-1) ) **2+ (Y (T) -Y (I- 1) ) **2) 

ARCL= S (NST) 

DO 60 1=1, NST 
60 SOL (I) = S (I) /ARCL 
C 

C CALCULATE P F FS , TJ E, M E , POPTZ , A N D VOVCR AT EACH STATION 
C 

C VELOCITY OVER CRITICAL VELOCITY GIVEN AS INPUT 
ISO DO 160 1=1, NST 

POPTZ ( I ) = (1. — GMM1/GMP1*V0VCR (I) **2) ** (GAM/GMM1) 

UE (I) = SORT (2. *G A M/GKM 1 *PT7/RHTZ* (1. -POPTZ (I ) ** (GMM1/GAM) 

*) ) 

TSE(I) = TTZ-UE (I) ** 2 / (2. *CP) 

AE (T) = SQRT (GAM*R*TSE (I) ) 

M E ( I) = UE(I)/AE(I) 

160 PRES ( I) = PCPTZ (I) *PTZ 
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C PRINT INITIAL CALCULATED PARAMETERS 
C 

170 WRITE ( 0 , 1000) 

WRITE (6 r 1010) PSZ.TSZ, HZ , AS 7., ATZ ,RHSZ, PHTZ, MUSE , MUTZ, NUSZ. NUT7,CP, 
1PR,TC,ARCL 

WRITE (6 , 102 C) (I, PRES (I) ,UE (I) , ME (I) , PCPTZ (I) .VOVCR (T) ,1=1, NST) 

C 

C 

C PRINT GEOMETRY PARAMETERS 

r* 

v- 

200 IF (KPRE.NE.1) GO TO 210 

WRITE (6, 103C) (I, X (I) , Y ( I ) ,S (I) ,XCM (I) # YCM (I) ,SOL(T) ,1=1, NST) 

C 

C CALCULATE OTHER NECFSSARY PARAMETERS AT EACH STATION 
C 

210 no 220 I -1 , N ST 

TEM 1 = 1. +GMK1/2.*MF (I) **2 
PHSW (T) = PRES (I) /R/TWAL (T) 

P HS F ( I) = PRES f I ) /R/TSE (I) 

HEADW f I ) - . 5*RHSW (I) *OE (I) **2 
HEADE (I) = . 5*PHSE (I) *UF (I) **2 
SW (I) = TKAL (I) /TTZ-1. 

SUTHL (I) = SORT (TKAL (I) /TTZ) * (^TZ + TCON) / (TKAL (I) + TCON) 

NUW(I) = SUTBL (I) *NUTZ* (1 . +SK (I) ) **2*TEM 1** (OAM/GMM 1) 

R W ( I ) = UE (I) *S (I) /NUK (I) 

TAWL (I) = TSE (I) * (1. + PR** ( 1 . /2 . ) * (T f M 1 - 1 . ) ) 

TAWT (I) = TSE (I) * (1.+PR** ( 1 . / .7 . ) * (TEM1-1.) ) 

THAR (I) ■= .5* (TWAL (I) +TSE (I) ) +.22*PP** ( 1 . / 3 . ) * (TTZ-TSE (I) ) 

MUBAE (I) = MUTZ* SUTHL (I) *T3AP (T) /TTZ 

BB (I) = ME (II *ATZ/NUTZ* (TSE (I) /TTZ) ** (GMP1/ (2.+GMM1) ) 

A A (I) = BB (I) *TSE (I) /TBAR (I) * (MU BAP (I)/MI!T7) **. 268 

FF(T)=1.+.159P*HE(I)**2 + .60*SK(I)+.2101*SW(I)*HF(7)**2+.0im*ME(I) 
1**4+_P180*SW (I) * M E (I) * * 4 + . 1R25*SW (I) **2+.0735*SW (I) **2*ME (I ) **2 
2+. 007 3* SW (I) **2*KE (I) **4 
220 CONTINUE 
C 

C FINITE DIFFERENCE METHOD USED TO CALCULATE VELOCITY AND MACH NUMBER 
C GRADIENTS ALONG THE SURFACE 
C 

DUDS ( 1 ) = (UE (2) - UE { 1 ) ) / (S (2) - S(1)) 

D M D S ( 1 ) = (ME ( 2 ) - KE(1) ) / ( S ( 2) - S(1)) 

DO 230 1=2, NST 
IM = I - 1 

IF (I .EQ. NST) GO TO 230 
IP = I + 1 

DUDS (I) = (UE (TP) -UE (IM) ) / (S (IP) -S (IM) ) 

DM DS (I) = (MF (IP)-ME (IM) ) / (S (TP) -S (TM) ) 

230 CONTINUE 

DUDS (NST) = (UE (NST) - U E ( I M) ) / ( S ( NST) - S(IM)) 

D MDS (NST) = (M E (NST) - M F ( I M) ) / ( S ( NST) - S(TM)) 

240 DO 2 5 0 1 = 1, NST 

250 DMDL (I) = ARCL*DMDS(I) 

r* 

C PRINT CT H F R CALCULATED PARAMETERS 
C 

IF (KPPE. NE. 1) GO TO 200 

WRITE (6, 105C) ( I , A E ( I ) ,TSE(T) ,TKAL(T) ,TA'/L (T) , T AWT (I) r t p a p (T) , 

11=1, NST) 
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WRITE (6 , 106 0) (T,FW (T) ,SW(I) ,S!!THL fl) ,PHSR (I) ,PH.SF (T) , HE ADW (I) , 

1 H E A D E ( I ) , N II Vi (I) , MUEAR (I) ,1 = 1 , I! ST) 

260 IF (KGFAD. NE. 1) GO TO 270 
WRITE (6 , 1 070) 

WRITE (6, 1C80) (I,DtJDS(T) ,DMDS (T) ,DMD1 (I) ,1=1, l?ST) 

CHECK FCP IMPROPER IN^UT 

270 DO 280 1=2, NST 

IF (UF(I) . N E . 0 . ) GO TO 280 
F P ROR = .TRUE. 

WR ITE (6 , 1090) 

RETURN 

280 CONTINUE 
RETURN 

FORMAT STATEMENTS 

1C00 FORMAT (1H1///4X, 24HPHFLIMINARY CALCULATIONS///) 

1C10 FORMAT (5X , 1 CHPSZ = F 1 2 . 5/8 X , 1 0 UTS Z = F10.4/5X, 10HUZ 

1 FI 1. 5//5X, 10HASZ = F 1 1 . 4/5X , 1 OH ATZ = F 1 1 . 4//5X , 1 0 HP HSZ 

2= G15.7/5X, 10HRHTZ = G1 8. 7//SX, 1 0HMUSZ = G 1 5 . 7/5 X , 1 0 H MUT7. 

3 = G 15. 7//5X, 10HNHSZ = G 1 6 . 7/5 X , 1 OHNUTZ = G 1 5. 7//5X, 10HCP 

4 = F11.4/5X, 1CHPR = F9 . S/5 Y , 1 0HTC = G1 5. 7/Ex , 1 OUARCL 

5 = F9.4///) 

1020 FORM AT (/IX, 7HSTATION,7X, 4HPRFS, 1 3X, 2HIIF, 12X ,2HMF, 1 1 X,?HP0PTZ,4X, 5H 
1VCVCF/(2X,I3,5X,F12.5,3X,F12.5 ,4X,F10.6,4X,F10.6,4X,F10.6)) 

1 0 30 FORMAT (///1X,7HSTATION,7X, 1 HX , 1 2 X , 1 H Y , 1 2 X , 1 US , 1 2X , 3 H XOM , 9 X , 3H YO M , 
19X,3HSOL/(2X,I3,3X,F12.S,1X,F12.S,lX,F12.5,4X,F9.S f 3X,F9.5,3X,F9.5 
2) ) 

1080 FORMAT (///1X,7HSTATION,5X,2MAE, 10 X, 3 HTSE , 9X , 4HTW A L , 8X , 4 II T A W I. , 8 X , 4 H 
1TAWT,8X,4HTRAR/(2X,I3,4X,F9.3,5(4X,F8.3))) 

1 06 0 FORMAT (///I X, 7HSTATI0N, 1 1 X , 2 8 R VI , 6 X , 2H S W , 4 X, 5HSUTHI , 7 X , 4 H RH S W , 1 2 X , 4 
1HRHSE,8X,5HHEADW,4X,5HHFAf. E,9X,1HNUW,12X,5HMUBAP/(2X,I3,3X,F16.1,2 
2X, F4. 1, 1X,F7. 3, 2X,G14. 6, 2X,G14. 6, 1X,FP. 3, IX ,F8. 3,2X,G14. 6,2X,G14.6 
3) ) 

1 C70 FORMAT ( 1H1///21X, 17HSURFACF GRADIENTS///) 

1080 FORMAT ( 1X,7HSTATI0N, 13X, 4HDUDS, 1SX, 4'ir-MDS, 1 5X,4HDMDL/ (2X,T.3 ,4X., F18 
1.6, 1X,F18.6 , IX, FIR. 6) ) 

1090 FORMAT ( // /// 1 0 X , 8 3 HT H E P E IS A STAGNATION PC TNT AT A STATION OTHFR 
1THAN STATION 1. THIS IS MOT AIIOWFE) 

END 

SIP FTC LAMNA LIST, DECK 
SUBROUTINE LAMNAR 

CO M MCN/G AMP E/G M P 1 , G M M 1 
CO MMON/CTOPL/DUM (2) , X M A X , DU M2 ( S ) 

COMMON/C1/G AM, R, PT7,, TT7 , UPM ACH , NST , NVP, NTI’PP.,KPVM, KFM, KSMTFF , 

1 KSPI. N , KLF , K ATCH , CTHFT , DI AM ,TL AM , DTUPP ,TTUPP ,KPPE,KGPAD, KSDF, KLAM , 
2KMA IN, K PROF, X (100) , Y (10°) , PPF.S ( 100) , UF (IOC) , ME ( 100) , POPT7. ( 100) , 

3 VO VCR (100) , TKAL (100) (10 0) 

CO MMCM/C2/PSZ , TSZ , UZ , AS 7 , AT7 , PFiSZ , P1IT7 , M US7 , M UT Z , NU SZ , NUT 7. , C P , 

1 PR , TC , ARCI. 

COM MON/C 3/ XOM (100) , YOM (100) ,S (100) , SOI. (100) , AF ( 100) ,TSE ( 100) , 

IT AWL (100) , TART (100) , TEAR (100) , RW (100) ,SW (1^0) , SUTHI. ( 1 00) , 

2RHSK (100) , RHSF. (100) , HFl'PK (100) , HFA DE (100) ,i:uv (IOC) , MU n AP (100) f 

3 A A (100) , PB (100) , FF (100) , DUDS (100) , PUDS (1°0) , D M D I (100) 
C0.MM0N/C4/TUFT ( 10O) , DELSP ( 1 00) , D FI.T?. ( 1 OP) , frp (KOI , 

1FOPMT (10 0) ,FCPMTR (100) , FTM (100) , RTF T (1^0) ,CF (190) f 
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1TAUVJ (100) , NUSS (100) ,DTDY (10G) , HTBAN (100) , CRN (100) 

CO MMON/C5/SKAPL (100) , SH A FK (100) ,E,NS 
COMMON/C6/FTR AN, POP MS 
COMMON/C7/I NST, ITR AN , ISFP 
COMMGN/C9/EEROR,TPANS, SEPRN 

DIMENSION CCPIN (100) ,COFML (100) , SHEAR (100) , DTH (100) 

DIMENSION CCN (20) , CP CP (20) ,CDIF (20) ,CSFR (20) ,CCRN (20) ,CDTH (20) 
DIMENSION STAB (505) , CTAP1 (505) ,CTAF2 (505) 

BEAL MUSZ,NUSZ,MUTZ,NUTZ,ME,NUW,MUEAP,NIISS,NUPK,KBAR, INTI, I NT 2 
LOGICAL ERROR, TRANS, SEPFN 
EXTERNAL FII NOT , I NT 1 , I NT 2 

READ DATA FOR CORLN(I), RCPIT, DTFF, SHEAR, CRN, AND DTH CU»VF FITS 

DATA (CCN (I) ,1=1,6) /- .08178, .06670, -.03143, 

I. 00 87 3,. 01 6 57, -.01 0 52/ 

DATA (CRCR (I) ,1=1, 6) /5. 47071, 43. 6053,227. 198, 

1-2067.04,-27172. 7, 13691.2/ 

DATA (CDIF (I) ,1=1,6) /903 . 785 , 26 3 65 . 0 , 3. fi 5695E+5, 

II. 1 10 4 4E+6.-4. 5 3853E+7,-7.702 76E + 7/ 

DATA (CS HR (I) ,1 = 1,16) /. 2244 88,-1.91539,-9.894,-68. 13488, 

1-. 001512,-1. 4768,- 10. 529 25,- 152. 2 781,-. 0024 06, -.01562°, 

1-1. 45743,- 126. 23395, . 000752,. 005385,. 917838,-39. 40644/ 

D AT A (CC RN (I), T= 1, 1 6) /2. 02056,-1°. 72 1 1, -2 4. 04 95,-1400.002, 

1-. 050979, -10. 880 12, 6 2. 4419,-5081 .76,-. 0 14343,2. 27°845, 
1129.7008,-6257. 848, . 0270567,- 1 . 67705 1 , c 7 . u 3 97 255 2 . 2 66/ 

DATA (CDTH (I) ,1 = 1 , 16) /8. 0 28 29, -4. 3 09 78, 88. 8 24 4, 36. 43 36, 
12.71101,-7.42259,242.293,-16.293, -.16394,-7.61843,286. 9795, 

164. 11 186,-. 16758,-3.70289,1 30.8107, 1 1 1.3276/ 

INITIALIZE PARAMETERS 

INST = 0 
ITRAN = 0 
ISEP = 0 
CF(1)= 0. 

TAU»(1)= 0. 

NUSS (1) = 0. 

DTDY(1) = 0. 

HTR AN ( 1 ) = 0. 

CRN (1) = 0. 

RTRAN= 0. 

IF (CTHET .GT. 0.) PATCH = 1 

CHECK CONSISTENCY OF INITIAL VALUES 

IF (DLAM.GE. 0. . A ND . TL AM . GE . 0 . . A N D. PTUPB. C-E. 0 . . A ND. TTUR B. GE. 0. ) 

1GO TO 10 
ERROR = .TRUE. 

WRITE (6,1000) 

RETURN 

10 IF (NTURB.NE.1) GO TO 30 
ITRAN = 1 

IF (DTURB.GT. 0. . AND. TTURB. GT. 0. ) GO TO 20 
ERROR = .TRUE. 

WR ITE (6 , 1010) 

BFTUPN 
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20 IF (UE (1) .GT.O.) GO TO 240 
ERROR = .TRUE. 

WRITE (6 , 1020) 

RETURN 

BEGIN CALCULATION IN LAMINAR REGION - CHFCK FOR INITIAL VALUES 
CALCULATE INITIAL CORRELATION NUMEER 

30 IF (DLAM.EQ. 0. . AND.TLAM. EQ. 0. ) GO TO 70 
IF (UE (1) -GT.O.) GO TO 40 
ERROR = .TRUE. 

WRITE (6, 1030) 

RETURN 

40 IF (TIAM.EQ.O.) GO TO 50 
INITIAL MOMENTUM THICKNESS WAS GIVEN 
TEM 1 = 1. +GMM1/2. *ME ( 1) **2 

COR ML (1)= -ATZ*TL AM**2/NUTZ/SUTHL (1) /AFC L/T EMI** ( (3. -GAM) / 

1 (2. *GMM 1) ) 

CCRLN(I) = CORML (1) *DMDL (1) 

GO TO 90 

INITIAL DISPLACEMENT THICKNESS WAS GIVEN 
50 IF (ABS (DMDL (1) ) .GE. .0001) GO TO 60 
COPLN (1) = 0. 

TEM 1 = 1.+GMM1/2.*ME (1) ** 2 

FORM (1)= 2. 38411* ( 1 , + (2.79-1.7R*PR**.5)* ( (1. + SW (1) ) *TEM1-1 .) ) + (4.6 
1 5 * P R * * (1./3.) -3. 65* PR*. 5) *PF**.5* (TEM1-1.) 

THET ( 1 ) - DLAM/FORM (1 ) 

CORML (1) = - ATZ*THET (1) * * 2/N UT Z/SUTH L (1) /AFCL/TEM1** ( ( 3 . -G A M) / (2 . * 
1GMM1) ) 

GO TO 90 

60 IF (DMDL (1) .GT.O.) CALL ROOT B (- 1 ., 0 DLA M , FU NCT 5E-5, COPLN ( 1) , SL) 
IF (DMDL (1) , IT. 0. ) CALL ROOTB (0 . , . 2 , DL A M , FUNCT, . 5E-5,CORLN (1) ,SL) 
CORML(I) - CORLN (1) /DMDL (1) 

GO TO 90 

MO INITIAL LAMINAR VALUES GIVEN 
CALCULATE INITIAL CORRELATION NUMEER 

SHARP LEADING EDGE 

70 IF (KLE. NE. 1 . AND. ABS (DMDL (1) ). GE. . 0001) GO TO RO 
CORLN (1) = 0. 

CORML (1) = 0. 

GO TO 90 

STAGNATION POINT 

BO CALL CURVFT (CCN, CORLN (1) ,SW (1) ,0,5,0) 

CORML ( 1 ) = CORLN (1)/DMDL (1) 

IF (CORML (1) . LT. 0.) GO TO 90 
EFROR= .TRUE. 

WRITE (6, 1040) 

RETURN 

SOIVE LAMINAR DIFFERENTIAL EQUATION 

CALCULATE CORRELATION NUMBERS ALONG THE SURFACE 

90 TEM1 - 1. +GMM1/2. *ME ( 1) **2 
TEM 2 = (3.*GAM-1.)/(2.*GMM1) 

DEL= 0. 002* ARCL 
SS= -DEL 
NTAB=1 
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CTAB1 (1) = CCBLN (1) 

CTAB2 (1) = CCRML (1) 

STAB (1) = 0. 

100 SS= SS+DEL 
SSDEL = SS+EEL 

CALL LGRNGE (S,SW,NST,SS,ANS1) 

CALL LGRNGE (S, ME, NST, SS, ANS2) 

CALL LGRNGE (S , M E , NST , -S S C EL , A NS 3 ) 

CALL LGRNGE ( S , DM DL , N ST, SSDEL , ANS4) 

A 1= 0.43631-0.00367+ ANS 1 +0. 0 0 4H 1*A NS 1**2 + 0.00651* ANSI** 3 
A 2= 5. 4 3 220 + 2. 254 no * A N Si — 0.06 6 7 2* AN Si** 2-0. 206 37* A NS1**3 
A 3= 4.51903-10. U9775 + ANS 1- 1 2. 7 1 7 32* ANS 1 **2- 2. 95 270* ANS 1 **3 
A 4- 19. 018 31+62.765 97* A NS1+ 1 15. 009 86* ANS 1** 2 + 62. 53 1 13* ANS 1**3 
A= A 1- A 3*CT AB1 (NTAB) **2-2.*A4*CTABl (NTAB) **3 
B= A2 + 2.*A3*CTAB1 (NTAB) +3.*A4*CTAB1 (NTAB) **2 
C 

C FOR SH - 0.0 

TF { CTAP1 (NTAB) .GE.-. 1) GO TO 101 
A=. 3953 
B= 4. 7 39 

101 K1 = 0 

SOLI = SS/ A PCI. 

SOL2 = SSDEL/ A PC L 

TEN 3 = SIMPS1 (SOL 1 , SOL 2 , INT 1 , K 1 ) 

IF (TEM3.EQ. C. .OR. K1 . EQ. 0) GO TC 110 
ER?OR= .TRUE. 

WRITE (6,1050) 

RETURN 

110 IF (NTAB.GT. 1) TEN4 = A NS 2 ** (- B) *T E M 1 **T F M 2 
TEM 1 ■= 1. +GMM1/2. *ANS3**2 
TE”I5 = ANS3** (-B) *TEM1**TEM2 
TEM6= -A*TEN5*TEM3 
IF (NTAB. EO. 1) T E M 7 = 0 . 

IF (NTAB.GT. 1) TEK7 = T E M 5 /T EK4* CT A E2 ( NT A B) 

NT AB= NTAE+1 
CT AB2 (NTAB) - TEM6+TEM7 
CT AB 1 (NTAB) = C? A B 2 ( NT A. B) * A N S 4 
STAB (NTAB) = SSDEL 
C 

C WHEN SW IS NOT EQUAL TO 0.0 , CURVE F TT RANGE ON CORLN 
C IS FRO* -. 32 TO . 16 
C 

IF (CT A El (NT A E) . GT. .50) GO TO 120 
IF (SS.LT.APCL) GO TO 100 
120 IF (KSDE.NE.1) GO TO 130 
WRITE (6, 1060) 

WRITE (6 , 1 07 C) (STAB (I) ,CTA31 (T) ,1=1 , NTAB) 

CALCULATE LAMINAR BOUNDARY LAYER PARAMETERS AT EACH STATION 

130 TF (KLAM.NE.1) GO TO 140 
WRITE (6 , 1080) 

140 1= 0 
150 1= 1+1 

IE (I.EO.NTURE) I TR A N=- 1 
IF (S (I) . LE. STAB (NTAB) ) GO TO 160 
WRITE (6,1090) 

160 IF (KLE .EC. 1 .AND. I .EO. 1) GO TO 151 
CALL LGRNGF (STAB, CT A B 1,NTAB,S (I) , CORLN (T) ) 

CALL LGRNGE (STAB, CTA.B2, NTAB, S (I) , CO R M I. (I) ) 
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C OETAIN SHEAR, CPN, AND DTH FPOM CURVE FITS VS COELN AND SN 
. 151 CALL CUR V FT (CSHP, SHEAR (I) , COR I N (I) , S W (T) ,3,3) 

CALL CURVFT (CCPN, CRN (T) , CORLN ( I) ,SW (I) ,3,3) 

CALL CURVFT (CUT H, DTH (I) , COP. IN (T) ,SW (I) ,3,3) 

C 

C FOE SW = 0.0 

IF (COFLN ( I) . OF. -.1) GO TO 161 
SHEAR (T) = -1. 2222+CCFLN (I) +.26 

CRN ( I ) = -58. 824 *CORLN (I) - . 676 5 

DTH (I) = -22. 222*COFLN (I) +7. 1 112 
C CALCULATE OTHFP LAMINAE BOUNDARY LAYER PARAMETERS 
161 T EM 1 = 1. + GMM 1/2. *MF, ( T) **2 

THET (I) = SQPT (-CORML (I) *NUTZ*S!JTHL (I) * A PCI,/ AT7,*TEM 1** ( (3. -GAM) / 

1 (2. *GMM 1) ) ) 

FORM (I) = (-1.1 1 3 8*CORLN ( I) + 2. 38 U 1 1 ) * ( 1 . + (2 . 70- 1 . 78 * PP * *. 5) * { (1 . + 
1SV (I) ) *TEM1-1 .) ) + (4. 65*PP** (1./3. ) - 3. 65*PP**. 5) 5* (TEM1-1.) 

DELSR (I) = THFT (I) *FORM f I) 

RTH (I) = LIE (I) *THET (T) /NUB ( I ) 

FORM (I) = (FORM (I)— SORT (PR) * (TEN 1-1.) )/( (1.+SK (T)) *TE M 1 ) 

FOR MT P (I ) = FOP MI (I) * (1.+SW (I) ) 

DELTA (I) = THET (I) * (DTH (I) + (TFM 1- 1 . ) * (FOPMTR (I) + 1 . ) ) 

SHAPL (I) = DELTA (I) ** 2/N U W ( J) * DUDS ( I ) 

IF (I.EQ.1) GO TO 180 

CFRK= 2.*SHEAP (t) *SORT (-SOL (T) /ME (I) /COD Ml (n ) 

C F (I) = CFRW/SORT (RW (I) ) 

T AUW ( T) = CF (I) *HEADW (I) 

NUFV= CFFW*PR**. 3/C«N (I) 

NUSS (I) = N tJ F W * S 0 F T (PW (I) ) 

DTDY (I) = NUSS (I) * (T AWL (T) -THAI. (I) ) /S (I) 

HTR AN (I) = TC*DTD Y (I) 

IF (TAUW (I) . GT.O. ) GO TO 180 
IF (KATCH.NF.O) GO TO 170 
ISEP= I 
SEPPN= .TRUE. 

RETURN 

170 ITRAN= -2 
GO TO 270 

180 IF (I. EC. 1. AND. UE (1) . EQ. 0.) GO TO 1°Q 

SHAPK (I) = NUTZ* RTII (I) **2*SUTUI (I) **2*(1.+SW (T) ) **U/ATZ. / ME (T) **2/ 
IFF (I)/AFCL*EMDL (I) * T E M 1 * * (1./GMM1) 

CO TO 200 
180 SHAPK (I) - 0.07 

200 FTHI (I) = RTH(I)*SnT!IL(I)*(1.+SW(I))**?/FF(I) /SORT (T'F.M 1 ) 

C 

C CALCULATE RCRIT TO CHECK FOP INSTABILITY AND TRANSITION 
C 

CALL CURVFT (CPCR , PCR IT , S HA PK (T) ,0,6,0) 

IF (SHAPK(I) .GT. .07) PCPIT - 8.3163 
RCRTT= FXP (RCRIT) 

IF (INST. NE. 0) GO TO 210 
C 

C CHECK FOR INSTABILITY 
C 

IF (PTHI (I) . LT. RCRIT) GO TO 270 
RINS= RTH 1(1) 

I N ST= I 
GO TO 270 
C 

C CHFCK FOR TRANSITION 
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210 K 1 = 0 
N 53 = I 

TE M = SIMPS 1 (SOL (INST) , SOL (I) , I NT 2, K 1) 

IF fTFM.Fq.C..OT>.K1.EO.O) 00 TO 220 
FPR0R= -TRUF. 

WRITE (6, HOC) 

RETURN 

220 KBAP= TEN/ (SCL (I ) -SOL (INST) ) 

CALL CURVFT (CCI F, DTFF, KPAR, 0,5,0) 

IF (FBAR . GT. -03) DIFF = 44f>00.*K EAF+700.0 
R T R A N = RINS+E1FF 
IF (RTHI (I) . LT. RTPAN) GO TO 270 
IF (I.LT.NTURE) GO TO 270 
ITFAN= -1 
GO TO 270 
230 IT P AN= I 

COMPOTE INITIAL VALUES FOR TUP ROLF NT SOLUTION 
2 UO TRANS= -TRUE. 

IF (DTURB. E0.0. , AND. TTUPE. F0-0.) GO TO 260 
IF (DTURB. GT.O.. AND. TT UP3.GT.0.) GO TO 250 
ERROR = .TRUE. 

WRITE (6,1110) 

RETURN 

250 THET (ITR AN) = TTURB 

FORM (IT RAN) = FTURP/TTURB 

TEM 1 = 1.+GMMl/2.*ME (ITRAN) * * 2 

FOP KT (ITP AN) = (FC R M (ITRAN) - PR** ( 1 . / 3. ) * (TEM1- 1 .) ) / ( (1.+SW (ITR AN) ) 
1*TFM 1 ) 

260 IF (CTitET.GT. 0. . AND. DTURB. EQ.O. . AND. TTURB. FQ.O. ) THET (IT RAN) = 
1CTHET*THET (ITRAN) 

THETTP = THET (ITRAN) * (TSF (ITRAN) /TTZ) ** (GMP 1/ (2. *GNM1) ) 

FTP AN = (ME (ITRAM) * A T 7*T H FTT F/N UT Z) * * 1 . 26 8 
IF (RTRAN .LE. 0.) GO TO 265 

FO R MS = FORMI (ITRAN) - 0 . 5 0 .3 <3Q-0. C6501*ALCG ( PT P A N ) +0.001272* (ALOG (PTF 
IAN ) ) **2 

IF (DTURB. GT. 0. .AND. TTURR.GT.O.) FCRMS=FO R MI (ITRAN) 

RETURN 

265 FORMS = 1.4 
RETURN 

PRINT OUTPUT 

270 IF (RLAM.NE. 1) GO TO 2 n 0 

IF (INST. FQ. 0 .OR. INST- EG. I) W P TT E ( 6 , 1 1 2 0 ) I ,CORLN (T) , SHEA R ( I) , 

1 FT H ( T) , FORMTR (I) ,SHAPL (I) , RTH I (T) .SHAPE (I) , PCRJT 
IF (INST. NE. 0 -AND. INST. NS. T) W F IT F (6 , 1 1 3 0 ) I , CORE N ( I ) , SHE A F ( I ) , 
1DTH (I) .FORMTR (I) ,SHAPL (I) , RTHT (T) .MPA 0 , DTFF, RTRAN 
IF (ITRAN. FC. -2) WRITE (6 , 1 140) 

2«0 IF (ITRAN. EQ. -1. OR. IT»A". *0.-2) 00 TO ?3n 

IF (I.EO.NST) PFTURN 
GO TO 150 

FOFMAT STATEMENTS 

1000 FORMAT {/////, 1 OX , 601! A NEGATIVE INITIAL VALUE HAS n E F N GTVFN. THIS 
IIS NO’ 1 ’ A I LOW ED) 
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1010 FORMAT {/////, 10X, 7512 INITIAL VALUES WERF NOT GIVEN FOR THE THRB'JLFN 
IT BOUNDAFY LAYER AT STATTON 1) 

1020 FORMAT (/////, 10X, BOH INITIAL VALUES WERE GIVEN FOR THE TURBULENT BO 
1UNDAFY LAYER AT A STAGNATION POINT) 

1030 FORMAT (/////, 1 9 X , U4 H IN IT I A L VALUES OTHFR THAN ZERO WERE GIVFN FOP 
1THE LAMINAR BOUNDARY LAYER AT A STAGNATION POINT) 

1 040 FORMAT (/////, 10X, 106HFOR THIS INPUT DATA STATION 1 IS ASSUMED TO ?. 
IE A STAGNATION POINT, STNCE NO INITIAL THICKNESSES ARE GTVEN./ 
210X, 1 18HIN THIS CASE PRESSURE SHOULD DECREASE INITIALLY. EITHER G 
3IVE AN INITIAL VALUE FOR DLSPLACFMENT OR MOMENTUM THICKNESS,/ 

4 1 OX , 6 OHOR BEGIN WITH A SMCPT REGION OF FAVORABLE PRESSURE GRADIENT 

5 .) 

1050 FORMAT (/////, 10X, 37HEPRCR IN COMPUTING INTEGRAL FOP CORLN) 

1060 FORMAT ( 1 H1///7X, 50HLAMINAR DIFFERENTIAL EQUATION - SOLUTION FOR CO 
1RLN///5 (24H S CCRLN ) //) 

1070 FORMAT [ < 5 ( F 1 2. 5 , 2X , r 7 . 4 , 3X) ) ) 

I 0B0 FORMAT (1H1///1 X,59HLAMTNAR CALCULATION OF INSTABILITY AND TP.ANSITI 

ION LOCATION S/// 1 X , 7HST ATION , 2 X , 5HCOH LN , 5 X , 5 HSHEA R, 5 X, 3 HDTH , 6 X ,6 HFO 
2RMTR,4X,5HSHAPL,9X,4HRTHI,6X, r )HSHAPK,PX,5HRCRIT,6X,4HK3AP,10X,4HDI 
3FF,9X,5HFTP AN) 

1090 FORMAT (/////, 1 OX, 65HLAMINAR SOLUTION HAS PROCEEDED BEYOND THE RANG 
IE WHFRE IT IS VALID) 

1109 FORMAT (/////, 1 0 X , 36HEPROR IN COMPUTING INTEGRAL FOR KBAR) 

1110 FORMAT (/////, 10X, 54HIF INITIAL TURBULENT VALUES APE GIVEN, THEY BO 
1TH MUST EE NONZERO) 

1120 FORMAT (14, 1X,5F10.4, 1 X , F 1 2 . 1, IX , F 1 0 . 5 , IX, FI 2. 1) 

II 30 FORMAT (14, IX, 5F10.4, 1 X , F 1 2 . 1, 24X,F12.5, 1 X , F 1 2 . 1 , 1X,F12. 1) 

1140 FORMAT ( ///// ,10X,85HLAMINAR SEPARATION HAS OCCURRED. ASSUMED TO BE 
1 TRANSITION TO TURBULENT BOUNDARY LAYFR) 

END 

SIEFTC TUPBL LIST, DECK 
SUBROUTINE TUPBLN 


COMMON/GAMPM/GiMPI , GMM 1 

COMMON/C 1/GAM, P , PTZ , TTZ , 11 P M AC H , N ST , N V F , NT UR H , K P V M , K EM , KS MTH , 

1 K SP L N , K L F , K A TC H , CTH FT , D L A M , TI. A M , D' r U P B , TT U p B , K PR F , K GF A F , K S D E , KL AM , 
2KMAIN,KPPOF,X (100) , Y ( 100) ,PRFS(100) , UF ( 1 9 0) , ME ( 10 0) , ROPTZ (1 09) , 

3 VO VCR (100) , TKAL (100) , FT A (100) 

COMMON/C2/PSZ,TSZ,UZ,ASZ,ATZ,PHSZ, P H T Z,MUSZ,M!JTZ,NUSZ,NUTZ,CP, 

1PP ,TC, ARCL 

COM MON/C 3/ XC M (190) ,YOM (100) ,S (100) ,SOL (10 0) , AE (100) , TSF. (100) , 

IT AWL (100) ,TAWT (100) ,TBAR(100) ,°S (10 0) , SW (100) ,SUTHL (100) , 

2RFSW (100) , R H S E (100) , HFACW (100) , H FACE (100) , NUK (190) , MU BAR (190) , 

3A A ( 100) , FB ( 100) , FF (1 00) , DUDS ( 100) , UMDS ( 100) , nt-IDL (100) 

COM MON/C 4/THFT (100) , PELSR (100) , DELTA (10 0) , FORM (100) , 

IFOR MI (100) , FCPMTR (1 00) , PTH ( 100) , RT H I (100) ,CF ( 19 9) , 

1TAUW ( 100) , NUSS ( 10 0) , DTDY ( 1 0 0) , HTR AN ( 100) ,CP N (100) 
COMMON/C6/FTRAN,FOPMS 
COMMON/CR/TNST, ITPAN,TSFP 

COMMON/CP/XT AB (505) ,YTS31 (505) , Y T A B 2 (505) ^'’’AB 
COMMCV/C9/FR PCR, TRANS, SFPHN 

REAL M.USZ, NUS7 , MUTZ,NUTZ, ME , NU W , MU P A F , NUSS 
LOGICAL FRRCR, TRANS, SEPFN 


SOLVE TURBULENT BOUNDARY LAYER DIFFERENTIAL EQUATIONS 
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C USING RUNG A— KUTTA 
CALI RUNKUT 

IF (KSDE. NE. 1) GO TO 10 
WRITE (6,1000) 

WRITE (6, 1010) (XTAB(I) , YTAB1 (T) ,YTAB2(I) ,I=1,WTAB) 

10 DO 5 1=1 ,NTAE 

IF (YTAB2 (I) -LE. 2.S) GO TO 6 
WRITE (6,100) 

100 FORMAT (9H0 (TUPBLN) , 5X , 63HINCC MPRESSI ELE FO*M FACTOR GREATER THAN 
*2.8 - CASE TERMINATED) 

STOP 

5 CONTINUE 

CALCULATE TURBULENT BOUNDARY LAYFF PARAMETERS AT EACH STATION 
DO 30 I = ITR AN , NST 

IF (S (I) . LE- XTAB (NTAB) ) GO TO 20 
ISEP = 1-1 
SEPPN= -TRUE. 

RETURN 

20 TEM 1 = 1. +GMM 1/2- *ME (I) ** 2 

CALL LGRNGE ( XT A B , YT A B 1 , NT AB , S ( I ) ,F) 

THETTR= N UT Z * F * * . 7886/ME (I) /ATZ 

THET(I) = THETTR* (TTZ/TSE (I) ) ** (GMP 1 / (2-*OMM1) ) 

RTK (I ) = DF (I) *THET (I) /NOW (I) 

CALL LGRNGF (XTAB, YT A B2 , NT A B , S (I) ,FOPMT (I) ) 

FO R MT R ( I) = FORMI (T) * (1. +SW (I) ) 

FORM (I)= FORMTR (I) *TEM1 + PR** (1./3. ) * (TFM1-1. ) 

DELSR (I) = THFT (I) *FORM (I) 

PQW ER = 2.0/ (FORMI (I) - 1- 0) 

IF (FORMI (I) .IT. 1.02) POWEP=10P. 

DELTA (I) = (1. +POWFR) *DFLSR (I) 

C F ( I) = 0. 246*EXP (— 1-561*FORMI (I) ) * (U E ( I) *7H?T (I) /NUTZ/ (TFMl** (1-/( 
1G AM- 1.) )))**(-. 26 8) *TS F (I) /TBAP (I)* (MUBAR (I) /MUTZ) ** (. 26 8) 

T AUW (I ) = CF (I) *HFADE (I) 

IF (I. EO. 1) GO TO 30 

HTRAN (I}= CF (I) /2./P 1 ?** (?./3.) *FHSF (I) *US (I) *CP* (7 A WT (I ) -T W A L ( I) ) 
DTD Y (I) = HTRAN (I) /TC 

NIJSS (I) = S (I) * DT DY (T) / (TAWT (I) -TSAI (I) ) 

CR N { I) = CF (I) *RW (I) /NUSS (T) 

30 CONTINUE 
RETURN 

1000 FORMAT (1 H1///5X, 62HTUR3ULENT DIFFERENTIAL FQUATTONS - SOLUTION FOR 
1 F AND FORMI///4 (31K S F FORMI ) //) 

1010 FORMAT ( (4 (F10.5, 2X,F8. 1,2X,F7. 4,2X) ) ) 

END 

$ I E FTC PROF I LIST, DECK 
SUBROUTINE PRCFIL 

CCMMCN/CCDE/KCDE, DELSRL, DELSFU 
COMMON/CTOBL/DUM (2),XMAX,DUM2(5) 

COM MCN/LNK 2/DU M 3 (3) ,CSTAR,DUM3S(2) , GST A PI , DUM 4 ( 3) 

COMMON/C1/GAM, D , PT Z- , TTZ , UPM ACH , NST, NVP,NTURB,KPVM,KF V ,KSMTH, 

1 KSPLN , KLE , K A TCH , CTHET , DL AM , m I, A M , DTUE B ,TTU PR , K n FE, KGR A D, KS PE , XI. AM , 
2KMATN, KPPOF,X (100) , Y (100) , PRES ( 100) ,UE ( 1 0 0) , ^ (100) , POPTZ (1 00) , 

3 VO VCR (100) , TWA L (100) , ETA (100) 

COM MON/C 3/XOM (100) ,YOM (100) ,S (100), SOL (ICO) , AF(IPO) , TGE(IOO) , 

IT AWL (10 0) , TAWT (100) , TRAP (100) ,RW (100) ,Svt (10 0) ,SUTHT, (10Q) , 

2PHSW (100) ,RHSE (100) , HE ADW (100) ,HFADE (100) ,NUW (100) , MUBAR (100) , 
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3AA(100) , BB(IOO) , FF(IOO) , DUDS (100) , DMDS(IOO) ,DMDL(100) 
C0MM0N/C4/THET (100) , DELSR (100) , DELTA (100) , FORM (100) , 

1FORMI (100) , FORMTR (100) , RTH (100) ,RTHI (100) ,CF (100) r 
1TAUW (100) ,NUSS (100) ,DTDY (100) ,HTRAN (100) ,CRN (100) 

COMMON/C5/SHAPL (100) ,SHAPK(100) , B , NS 
COMMON/C7/INST, ITRAN ,ISEP 
REAL HE , NUSS 

PRINT LOCATIONS OF INSTABILITY, TRANSITION, AND SEPARATION 

IF (KMAIN.NE. 1) GO TO 60 
WRITE (6, 1000) 

IF (INST.EQ. 0) GO TO 10 
WRITE (6, 1010) INST 
GO TO 20 
10 WRITE (6,1020) 

20 IF (ITRAN.LE.1) GO TO 30 
WRITE (6, 1030) 1TFAN 
GO TO 40 
30 WRITE (6, 1040) 

40 IF (ISEP. EQ. 0) GO TO 50 
WRITE (6, 1050) ISEP 
GO TO 60 
50 WRITE (6, 1060) 

PRINT LOCATIONS OF LAMINAR AND TURBULENT BOUNDARY LAYERS 
60 IEND = ITR AN— 1 

IF (IEND. EQ.-1. OR. IEND. EQ.O) IEND=ISEP 
IF (IEND. EQ- 0) IEND=NST 
IF (KM AIN. NE. 1 ) GO TO 70 
IF (ITRAN.EQ. 1) WRITE (6 , 1070) 

IF (ITR AN. NE. 1) WRITE (6 , 1 080) IEND 
IF (ITRAN.EC-0) WRITE (6,1 090) 

IF (ITRAN.EQ. 1) WRITE (6 , 1 1 00) ITR AN , TEND 
70 IF (ITRAN.LE.1) GO TO 80 
TEND = ISEP 

IF (IEND. EQ. 0) IEND=NST 
IF (KMAIN.NE.1) GO TO 90 
WRITE (6, 1 100) ITRAN, IEND 

80 IF (KODE .EQ. 0) DELSRU = D ELS R (I EN D) /X M A X 
IF (KODE .EQ. 1) DEI.SRL = DELSR (TEND) /XMAX 

WRITE (6,2000) 

2000 FORMAT (1 HI , 35X,29HDTMENSIONED BLADE COORDINATES , 36X, 24HANGLES OF 
♦ROTATION (DEG) //37X, 1 1H UNCORRECT ED, 1 4X ,9HCORRECTED , 15 X, 1 OHTRANSL AT 
$ED//1X,7HSTATI0N, 12X, 1 H X , 23X , 1 HY , 2 1 X , 5H YCOF R , 18X,5HYTR AN) 

SIGN = +1. 

IF (KODE .EQ.. 0) SIGN = -1. 

GST AR = GST ARI*XMAX/CSTAR 
DO 170 1=1, IEND 

YCORR = Y (I) ♦ SIGN*AES (DELSR (I) /COS (ETA (I) ) ) 

YTRAN = YCORR - GSTAR 
ETAD = 57. 2957796*ETA (I) 

WRITE (6,2001) I,X(I) ,Y(I) , YCORR , YTRAN, ETAD; 

2001 FORMAT (2X , 13 , 9X , 4 ( F9. 5 , 1 5X) , F9 . 4) 

170 CONTINUE 
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C PRINT CALCULATED BOUNDARY LAYER PARAMETERS 
C 

IF (KMAIN -NE. 1) GO TO 90 
WRITE (6, 1110) 

WRITE (6, 1120) (I, X (I) ,S (I) , BFLSR (I) ,T BET (I) , DELTA (I) , FORM (I) r 

IFORMI(I) , 1 = 1 , LEND) 

WRITE (6,1130) 

WRITE (6, 1 14C) (I,CF (I) ,TAUW (I) , RTH (I) , DTDY (I) , NUSS (I) , HTRAN (I) , 
1CPN (I) ,1=1, IEND) 

C 

C COMPUTE BOUNDS ON VELOCITY PROFILES 
C 

90 IF (KERCF.NE.1) RETURN 
WRITE (6, 1150) 

IF (ITRAN. NE. 0) GO TO 100 
IL 1= 2 
IL2= IEND 
IT 1 = 0 
IT2= 0 
GO TO 110 
100 IL 1= 2 

IL 2= ITRAN— 1 
IT 1 = ITRAN 
IT 2= IEND 

IF (IT1.EQ.1) IT 1 = 2 

C 

C CAICULATE AND PRINT LAMINAR BOUNDARY LAYFP VELOCITY PROFILES 
C 

110 N VP 1 = NVP+1 

IF (IL2.LT.IL1) GO TO 140 
DO 130 I=IL 1 ,IL2 
WRITE (6, 1160) I 
AAA= 2. +SHAPL (I) /6. 

BBB= — . 5*SH APL (I) 

CCC= -2.+.5*SHAPL (I) 

DDD= 1,-SHAPL (I) /6. 

DEL= DELTA (I) /FLOAT (NVP) 

YP= -DEL 

DO 120 J= 1 , N VP 1 

YP= YP+DEL 

ETAA = YP/DELTA (I) 

Y X M A X= Y P/X (NST) 

UUE = ( ( (DDD*ET AA + CCC) *ETAA+BEB) *ETAA + AAA) *ETAA 
U= UU E* UE (I) 

120 WRITE (6,1180) ET A A , YP , Y XM A X , U , UUE 
130 CONTINUE 
C 

C CALCULATE AND PRINT TURBULENT BOUNDARY LAYER VELOCITY PROFILES 
C 

140 IF(ITI.EQ.O) RETURN 
DO 160 I=IT 1 , IT2 
POWER= DELTA (I) /DELSP(I) -1. 

WRITE (6, 1 170) I, POWER 
DEL= DELTA (I) /FLOAT (NVP) 

YP= -DEL 

DO 150 J= 1 , N VP 1 

YP= YP+DEL 

ETAA = YP/DELTA (I) 
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YXMAX = YP/XMAX 

UUE = ETA A** (1. /POWER) 

U = UUE* U E (I) 

150 WRITE (6,1180) ETA A , YP , Y XM A X, CJ , UUE 
160 CONTINUE 
RETURN 

FOEHAT STATEMENTS 

1000 FORMAT (1H1///1X, 36HPRINCIPAI BOUNDARY LAYER INFORMATION///) 

1010 FORMAT (/10X,31HINSTABILITY OCCURS AT STATION ,13) 

1020 FORMAT (/10X, 26HINSTABILITY DOES NOT OCCUR) 

1030 FORMAT (/ 10X , 30HTR AN SITION OCCURS AT STATION ,l3) 

1C40 FORMAT (/10X,25HTRANSITION DOES NOT OCCUR) 

1050 FORMAT (/I OX , 30 HSEPAR ATION OCCURS AT STATION ,13) 

1060 FORMAT (/10X,25HSEPAEAT.ION DOES NOT OCCUR) 

1070 FORMAT (/10X, 37HLAMINAR BOUNDARY LAYER DOES NOT OCCUR) 

1080 FORMAT (/10X, 42HL AMlNAR BOUNDARY LAYER - STATIONS 1 TO ,13) 

1090 FORMAT (/1 0 X , 3 9HTUR BULENT BOUNDARY LAYER DOES NOT OCCUR///) 

1100 FORMAT (/10X, 35HTURBULENT BOUNDARY LAYER - STATIONS, 2X, 

113, 6H TO ,13///) 

1110 FORMAT (/1X,7HSTATI0N,RX, 1HX, 12X ,1HS, 12X,5HDELSF, 10X,4HTHET, 11X, 
15HDELTA,11X,«HF0FM,10X,5HFORMI) 

1120 F0PMAT(2X,I3,3X,2F13.6,F14.6,1X,F14.6,1X,F14.6,1X,2F14.U) 

1130 FORMAT (/// 1 X ,7HST ATI ON , 6 X , 2 HCF , 1 3X,4HTAUW, 1 IX, 3HBTH, 14X, 4HDTDY, 
113X,4HNUSS, 10X,5HHTPAN, 12X,3HCRN) 

1140 FORMAT (I5,F14, 5, 2X,F14.5,1X,F12, 1, 5X,F14.2,2X,F14.2,1X, 
1F14-4,2X,F13.3) 

1150 FORMAT (1H1///1X, 17HVELOCITY PROFILES///) 

1160 FORMAT (/I X , 7HST ATION , 1 X , 15 , 2 X , 7 HPECFILE/3 X , 7HY/PELT A, 9X , 

11 HY, 12X, 6HY/XMAX, 10X, 1HU, 12X,4HU/UE) 

1 170 FORMAT (/ 1 X , 7BST ATION , 1 X , 15 , 2 X , 7H PROFILE , 28 X , 2HN= , IX, F6, 2/3X,7HY/P 
1 ELT A , 9X , 1 RY , 12X, 6HY/XMAX, 10X, 1HU, 1 2X ,4HU/UE) 

1 180 FORMAT(1X,F8.4,2X,2G15.fi,2X,F9.2,6X,F8.4) 

END 

IEFTC RUNKT LIST, DECK 
SUBROUTINE PUNKUT 

RUNKUT SOLVES SIMULTANEOUS FIRST CFDER INITIAL VALUE 

ORDINARY DIFFERENTIAL EQUATIONS 

COMMON/C 1 /GAM, R,PTZ,TTZ, U PM ACH, N ST , N V P , NTU P B , KP V M , K FS , KS MT H , 
1KSPLN,KLE, HATCH, CTHET , DLAM ,TLAM ,DTUFR,TTURB,KPRE,KGRAD, KSDE, KLAM , 
2KMAIN,KPROF,X (100) , Y (100) , PRES (100) ,UE (100) , ME ( 100 ) , POPTZ ( 1 00) , 

3 VO VCR (100) , TW AL (100) , ETA (100) 

COMMON/C 3/XOM (100) , YOM ( 100) ,S (100) , SOL ( 100) , AE ( 1 00) ,TSE (100) , 
ITAWL(IOO) , T A WT (100) ,T3AP(100) , RW (100) ,SW (100) ,SUTHL(100) , 

2RHSW (100) , RHSE (100) , HEADW ( 100) , HEADE (100) ,NUW (100) , MUBAR (100) , 

3AA (100) ,BB (100) ,FF (100) , DUDS (100), DM DS (100) , DMDL (100) 

CO MMON/C6/FT RAN, FORMS 
C0MM0N/C7/INST,ITRAN,ISEP 

COMMON/C8/XT AB (505) , YTAB1 (505) , YTAB2 (505) , N T A B 
DIMENSION YY (2) ,PY(2) ,YINC(2) ,DCT(2) ,RUK(2,4) 

DOUBLE PRECISION XX , R X , Y Y , R Y , RUK , DEL , DOT , 

1TEM1 ,TEM2,TEM3,TEM4 ,TEM5 ,TEM6 

Real me, nuk, mubar 

SET DEL SPACING AND STORE INITIAL VALUES 
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DEL= 0.002*S (NST) 

10 YY (1) =FTRAN 
YY (2) = FORMS 
XX= S (ITRAN) 

N V=2 

NTAB = 1 
YTAB1 (1) = YY (1) 

YTAB2 (1) = YY (2) 

XT AB ( 1) = XX 
C 

C SOLVE FOR Y Y ( 1 ) AND YY (2) AT NEXT XX INCREMENT 
C 

C SAVE PREVIOUS YY{1) AND YY (2) 

20 DO 30 J= 1 , N V 
30 RY (J) = YY (J) 

PX= XX 

CALCULATE NEW YY(1) AND YY(2) 

DO 90 L=1,4 

C PUT DIFFERENTIAL EQUATIONS IN THE FCHM OF 
C FIRST DERIVATIVE = REMAINDER OF EQUATION 
CALL LGRNGE (S, HE, NST, XX, ANSI) 

CALL LGRNGE (S, SW, NST, XX, ANS2) 

CALL LGRNGE (S,AA, NST, XX, ANS3) 

CALL LGRNGE (S, BB,NST, XX, A N S4) 

CALL LGRNGE (S,DMDS, NST, XX, ANS5) 

CALL LGRNGE (S, TEAR, NST, XX, ANS6) 

TEM 1 - 1.+ (1.+ANS2) *YY (2) 

TEM 2= . 12 3* EXP (- 1.561*YY (2) ) *ANS3 

DOT (1) = 1,268* (— YY ( 1 ) /A NS 1 * A NS5*TEM 1 +T EM2) 

TEM 3= YY (2) * (YY (2) +1.) **2* (YY (2)-1.) 

TEM4= 1.+ANS2* (YY (2) *YY (2) +4. *YY (2) -1.) / ( (YY (2) + 1.) * (YY (2) +3.) ) 
TEM 5= (YY (2) *YY{2)-1.) *YY (2)/YY (1)* (. 12 3* EXP (-1. 561*YY(2) ) *ANS3) 
TEM6= (YY (2) *YY (2) - 1.) / YY ( 1 ) ** (. 78 86 ) * (, 0 1 1 * f Y Y ( 2) + 1 . ) * ( YY ( 2) - 1. ) 
1**2/YY (2) **2*TTZ/ANS6) *ANS4 
DOT (2) = — ANS5*.5/ANS1*TEM3*TEM4+TE M5-TEM6 
C APPLY THE PUNGA-KUTTA SCHEMF 
DO 40 J= 1 , N V 
40 RUK (J, L) = BEL*DOT (J) 

GO TO (50,80,70,90) , L 
50 DO 60 J=1,NV 
60 YY (J) = RY (J) + RUK (J,L) /2. 

XX= RX+DEL/2. 

GO TO 90 
70 DO 80 J= 1 , N V 
80 YY (J) = RY (J) +PUK (J,L) 

XX= PX+DEL 
90 CONTINUE 

C INCREMENT THE DEPENDENT VARIABLES TO OBTAIN NEW YY(1) AND YY(2) 

DO 100 J= 1 , N V 

YINC(J) = (RUK (J, 1) +2. *RUK (J, 2) +2.*P0K (J, 3) +FUK (J,4) ) /6. 

100 Y Y ( J ) = RY (J) +YINC (J) 

C 

C STORE NEW COMPUTED VALUES IN A TABLE 
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NT AB = NTAB+1 
YTAB1 (NTAE) = Y Y ( 1) 

YTAE2 (NTAE) = YY (2) 

XT A E ( NT A E) = XX 

IF (YTAE2 (NTAB) .GT. 2-8) BFTUFN 
IF (XX. IT. S (NST) ) GC TO 20 
RETURN 
END 

$IEFTC SPLIN DECK 

SUBROUTINE SPLINE (X , Y, N , DYDX, D2YDX2) 

SPLINE FITS A SPLINE CURVE TO X AND Y 

AND CALCULATES FIRST AND SECOND DERIVATIVES AT THE SPLINE POINTS 
END POINT SECOND DERIVATIVES EQUAL THOSE AT ADJACENT POINTS 

DIMENSION X (N) ,Y (N) ,DYDX(N) , D2YDX2 (N) 

DIMENSION G (100) ,H (100) 

G ( 1 ) — "I- 
H(1)= 0. 

N1 = N— 1 

IF (N1.LT.2) GO TO 20 
DO 10 1=2, N1 
A- (X (I)-X(I-I) )/6. 

B= (X (1 + 1) -X (I) ) /6. 

C= 2.* (A+B)-A*G (1-1) 

D= (Y (1+1) -Y (I) ) / (X (1 + 1) -X (T) ) - (Y (I) -Y (1-1) ) /(X (I) -X (1-1) ) 

G ( I) = B/C 

10 H (I) = (D-A*H (1-1) ) /C 
20 D2 YDX2 (N) = H (N 1 ) / ( 1 . +G (N 1 ) ) 

DO 30 1=2, N 
K= N+1-I 

30 D2YDX2 (K) = H (K) -G (K) *D2YDX2 (K + 1) 

DYDX (1) = (X (1) - X (2) ) /6. * (2.+D2YDX2 (1) +D2YDX2 (2) ) + (Y (2) -Y (1) ) / (X (2) 
1-X (1) ) 

DO 40 1=2, N 

40 DYDX (I) = (X (I)-X (1-1) ) /6.* (2.*D2YDX2 (I) + D2YDX2 (1-1) ) + (Y (I) - Y (1-1) ) 
1/ (X (I) — X (1-1) ) 

RETURN 
END 

SIEFTC LGRNGV DECK 

SUBROUTINE LGRNGE (X , Y , N , ARG , AN S) 

C LINEAR INTERPOLATION (REPLACES 4-POINT LAGRANGE) 

DIMENSION X (N) , Y (N) 

IF (ARG -GE. 0.80*X(1)) GO TO 2 
WRITE (6,100) ARG, X ( 1 ) , X (N) 

100 FORMAT (9H0 (LGRNGE) ,5X, 10HABSCISSA =,E13.5,27H IS OUT OF RANGE - 
* X(1) =,E13.5,5X,6HX (N) = , E 1 3. 5/ 15X, 27HEXTR APOLATED VALUE RETURNED 
*) 

2 NM = N - 1 
DO 3 1=2, NM 

IF (ARG -GT. X(I)) GO TO 3 
M = I 
GO TO 4 
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3 CONTINUE 

IF (ARG . LE. 1 . 20*X ( N) ) GO TO 5 
WRITE (6,100) ARG,X (1) , X (N) 

5 M = N 

4 A NS = Y (M) + (Y(M) - Y (M-1) )/ (X (M) - X(M-1))*(ARG - X (M) ) 
RETORN 

END 

$ I EFTC SIMP DECK 

FUNCTION SIMPS1 (X1,X2,FUNC, KSIG) 

DIMENSION V (200) ,H (200) ,A (200) ,B (200) ,C (200) ,P (200) ,E (200) 

LOGICAL SPILL 

DOUBLE PRECISION A NS , Q 

DATA TWC,THFEF, FOUR, THIRTY/2. 0,3. 0,4. 0,30.0/ 

DATA T , NM A X , NSIG/3. OE-5,200, 1/ 

C INITIALIZE FIRST ELEMENTS OF ARRAYS. 

V=X1 

H= (X2-V) /TWO 
A = FTJNC (V) 

B= FUNC (V + H) 

C=FUNC ( X 2 ) 

P=H* ( A + FOU P* E + C) 

E=P 

ANS=P 

N=1 

FR AC=T 

SPILL-. FALSE. 

10 TEST=ABS (FRAC+ANS) 

K = N 

DO 30 1=1, K 

C TEST MAGNITUDE OF 4TH ORDER ERROR IN THIS INTERVAL. 

IF (AES (E (I) ) .LF.TEST) GO TO 30 
IF (N.LT.NMAX) GO TO 20 
C GO TO FINISH IF STORAGE IS FILLED UP. 

SP ILL= . TRUE . 

KSIG=KSIG+NSIG 
GO TO 40 

C S U EDI V ID F INTERVAL AGAIN TO R-EBUCE 4TH ORDER ERROR. 

20 N= N + 1 

V (N)=V (I) + H (I) 

H (N)=H (I) /TWO 
A (N) = E (I) 

B (N) =FUNC (V (N) +H (N) ) 

C (N) =C (I) 

P (N)=H (N) * (A (N) + FOU R * B ( N) +C (R) ) 

H (I) = H (N) 

B (I) =FUNC (V (I) +H (I) ) 

C ( I) = A (N) 

Q=P (I) 

P (I)=H (I) * (A (I) + FOUR*B (I) + C (I) ) 

Q=P (I) +P (N) -0 
A N S=A N S + 0 
E (I)=0 
E (N) =0 
30 CONTINUE 

C TEST ALL INTERVALS AGAIN IF ANY WERE SUEriVICED THE LAST TIME. 
IF (N.GT.K) GO TO 10 
40 0=0.0 

DO 50 1=1, N 
50 Q=Q + E (I) 
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C TIGHTEN ERROR LIMIT IF TOTAL ACCUMULATED ERROR TOO LARGE. 

IF (ABS (Q/T) .LE.ABS (ANS) .OR. SPILL) GOTO 60 
FRAC=F RAC/TWO 
GO TO 10 

C FINISH CALCULATION. 

60 SIMPS 1 = (ANS+Q/THIRTY) /THREE 
RETURN 
END 

SIEFTC CURVF DECK 

SUBROUTINE CURVFT (COEF, ANS , X, Y , NX , NY) 

EVALUATE THE POLYNOMIAL FUNCTION, ANS=F(X,Y), USING COEFFICIENTS, COEF 

DIMENSION COEF (20) 

NX 1 = NX + 1 
NY 1 = NY+1 
ANS = COEF (1) 

IF (X.EQ. . 0. AND. Y. EQ. . 0) RETURN 
IF (Y.EQ..0) GO TO 10 
IF (X.EQ..0) GO TO 30 
GO TO 50 
10 DO 20 1=2, NX1 
20 ANS = ANS + COEF (I) *X** (1-1) 

RETURN 

30 DO 40 1=2, NYl 
K = (1-1) *N X 1 + 1 
40 ANS = ANS+CCEF (K) *Y** (1-1) 

RETURN 
50 ANS = .0 

DO 60 1=1, NYl 
DO 60 J= 1 , N X 1 
K = (I-1)*NX1 + J 

60 ANS = ANS + COEF (K) *Y** (1-1) *X** (J-1) 

RETURN 
END 

SIEFTC AFTMIX LIST, DECK 

SUBROUTINE AFMIX ( AL P H 1 , DELS , DEL P, THETS, THPTP , TE, SP , X MFS 1 ) 

DIMENSION XXX (100) ,YYY (100) ,SSS (100) 

COMMON/GAMPM/GMP1 ,GMM 1 

COMMON/C1/GAM,R, PTZ, TTZ , UPMACH,NST,NVF,NTURB,KPVE,KEM,KSMTH, 
1KSPLN,KLE,KATCH, CTHET, DL AM , TL AM , DTURB ,TTURR , KPRE , KGR A D , KSDE , KL AM , 
2KMAIN,KPROF,X (100) ,Y (100), PRES (100), UE(100),ME(100) ,POPTZ (100) , 
3VOVCR (100) , TNAL (100) , ETA (100) 

COM MON/C 3/XOM (100) , YOM (100) ,S (100) ,SOL (100) , AE (100) ,TSE(100) , 

1TAWL (100) , TAUT (100) ,TBAR(100) , PH (100) ,SW(100) ,SUTHL (100) , 

2RHSH (100) , RHSE ( 1 00) ,HEADW(100) ,HEADF(100) .NUH(IOO) .MUBAR(IOO) , 

3AA (100) ,BB (100) , FF (100) , DUDS (100) ,DMDS (100) , DMDL (100) 

COMMON/C4/THET (100) ,DELSR (100) , DELTA (100) , FORM (100) , 

IFOR MI (100) , FORMTR ( 1 00) ,RTH(100) ,RTHI(100) ,CF(100) , 

1TAUH (100) ,NUSS (100) , DTD Y (100) ,HTRAN (100) , CRN (100) 

EQUIVALENCE (X,XXX) , (Y, YYY) , (S,SSS) , (NST ,N) 

REAL ME 

KODE = 0 
WRITE (6,5) 

5 FORMAT ( 1 H 1 , 20X , 22H AFTERMIXING PROPERTIES) 

XALPH1 = ALPH1 

ALPH1 = ALPH1*0. 017453 
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2 VVCP1 = SQRT ( (GMP1/2.*XMFS1**2)/(1.+GMMl/2. *XMFS 1**2) ) 

XX=SP*COS ( ft L PH 1 ) 

DELSRT= (DELS*DELP) /XX 
THET A = (THETS*THETP) /XX 
DTE = TE /XX 

APS 1 = GMM 1/GMP1 *VVCR1**2 
A = 1 .O-DELSRT— DTE-THETA 
A 1 = 1 . 0— DELS FT— DTE 
IF (A .LE. 0-0) GO TO 16 
IF (KODE .EQ. 0) WRITE (6,1) 

1 FORMAT (1H0, 10X, 39HROTOF WITH NO BOUNDARY LAYER CORRECTION) 

C = ( (1.-AFS1) *GMP1/ (2. *GAH) + (COS (ALPH1) ) ** 2 * 

1 A*VVCR1 ** 2 ) / (COS (ftLPHI ) *A1*VVCR1) 

D = VVCR1*SIN (ALPH1) *A/A1 

VXVCP2 = Gft K*C/GMP1 -SORT ( (G AM*C/GMP 1 ) **2 
1-1. +GMM1/GMP1*D**2) 

GO TO 4 
3 KODE=2 

WRITE (6,10) 

10 FORMAT (1H0, 10X, 19HSUPERSONIC SOLUTION) 

VXVCP2 = GAM*C/GMP1 +SQRT ( (GAM*C/GMF1) **2 
1-1. +GMMl/GMPl*D**2) 

4 DENCH - GMM 1/GMP1* (D**2+VXVCP2**2) 

IF (DENCH. GE. 1.0) RETURN 

DENSR2 = (1. — GMM1/GMP1* (D**2+ V X VCR 2* * 2) ) 

1** (1./GMM1) 

DENSR1 = 1. /( (1. + GMM 1/2. *XMFS 1**2 ) ** (1./GMM1) ) 

PR 2 = DENSP2**G AM 

PT2PT0'= (BENSR1*VVCR 1*COS (ftLPHI ) *A 1 ) / (DENSR2* VXVCR2) 

PT0P2 = 1.0/ (PT2PT0*PR2) 

EBAR2 ■= ( (PT2PT0) ** f-GMMl/G A M) - 1 . ) / (PT0P2** 

1 (GMM 1/GAM) — 1.) 

ETAN = 1.0-EEAR2 

VVCR2 = SQRT (D**2 + VXVCR2**2) 

XM 2 = SQRT ( ( (2./GMP1) *VVCR2**2) / (1 .— (GMM 1/GMP1) 

1 *V VCR 2**2) ) 

T2TT0 = 1.— GMM 1/GMP1 *VVCR2**2 
ALPH2 = ATAN (D/VXVCR2) 

XMX 1 = XMFS 1*COS (ALPH1) 

XMX2 = XM2*CCS (ALPH2) 
ftLPH2 = ALPH2*57. 2958 

WRITE (6,6) XMFS1 ,SP,TE ,XM2,VVCR1 ,XMX1 ,XMX2 

6 FORMAT (1H0,7HXMFS1 = , F6. 4 , 2 X , 9 HSPACI NG =,F7. 6,2X,4HTE =,F7.5,2X, 
15HXM2 =,F6. 4,2X, 8HV/VCR 1 = , F6. 3, 2 X , 6HXM X 1 = , F6. 3 , 2 X , 6HXMX2 =,F6.3) 

WRITE (6,7) XALPH1, ALPH2 , PT2 PTO , PTO P2, T2TT0 . VVCR2 , E BAR 2 , ETflN 

7 FORMAT (1H0,6HALPH1=,F7. 3,2X,6HALPH2=,F7. 3,2X,8HPT2/PT0=,F7. 4,2X, 
17HPT0/P2=,F9.3,2X,7HT2/TT0=,F7. 4 , 2 X ,7HV/VCR 2= , F6 . 3 , 2 X , 6 H EB AR 2=, 
2F7.5,2X,6HETA-N=,F6. 4) 

IF (NODE .EQ. 2 ) RETURN 
IF ( KODE. EQ. 1 ) GO TO 3 
16 SP = SP + (DELS + DELP)/COS (ALPH1) 

WRITE (6,20) 

20 FORMAT ( 1 HO , 1 0 X , 36HROTOR WITH BOUNDARY LAYER CORRECTION) 

KODE = 1 
GO TO 2 
END 
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n n n 


$ I E FTC FUNC DECK 

SUBROUTINE F IJNCT (XX, FX , P FX , INF) 


COMMON/GAMPM/GHP1,GMK1 

COMMON/C 1/GAM, 3 , PTZ , TT Z , 11 PM ACH , N ST , NVF, NTURB,KPVM, KFM, KSHTH, 
1KSPLN,KLE,K ATCH.CTHET, DLAM , TL AM , DT U P E , TTURB , KFRE, KGRAP, KSDE , KLAM , 
2KMAIN,KPPOF,X (100) ,Y(1O0) , PRES (100) ,UE(100) , ME ( 1 00) , POPTZ ( 1 07) , 

3 VO VCR (100) ,TWAL (100) ,ETA (100) 

COM MON/C 2/P SZ,TSZ, IIZ , AS Z , ATZ , PH SZ , RHTZ, M USZ , M UTZ, NUSZ , N UT Z , CP , 

1 PR , TC , ARCL 

COMMON/C3/XCM (10 0) r YOM ( 1 00) , S ( 1 00) , SOL ( 1 00) , A E (1 00) , TSE (100) , 

IT AWL (100) ,TAWT ( 1 00) ,TBAP (100) , RW (100) ,SW (100) ,SUTHL (100) , 

2PHSW (100) ,R USE (100) , HEADW(IOO) , HEAr.E(IOO) , NOW (100) , MUE AP(IOO) , 

3AA (100) ,BB (100) ,PF (100) ,DUDS (100) , DMDS (100) ,DMDL (10 0) 

REAL MUSZ,NUSZ,MUTZ, NUT 7, MF,NUW ,MUEAR 

INF = 0 

B 1 = 1. +GBM1/2. *ME (1) **2 

B2~ 1. + (2. 79-1. 78*0?**. 5)* ( (1.+SW (1) ) * I! 1 - 1 . ) 

B3 = -NUTZ*SUTHL (1) * A RCL/ ATZ/EM CL ( 1) *R1 ** ( ( 3. -GAM) / (2. *GMK1 ) ) 

B4= -1. 1 1 38*B2 

85= 2. 38411*B2+(4.65*PR**(1./3.)-3.65*PR**.5)*PR**.5*(Bl-1.) 

FX = (B3*XX) **.5* (B4*XX+B5) 

IF (XX.EO.O.) GO TO 10 

DFX= .5*(B3*XX)**(-.5)*B3*(B4*XX+B5)+B4*(B3*XX)**.5 
RETURN 
10 INF = 1 

DFX = 1.E10 

RETURN 

END 

SIBFTC ROT B DECK 

SUBROUTINE POOTB ( A , B , Y , FUN CT ,T CLEF Y , X , DFX) 

ROOT FINDS A ROOT FOR (FUNCT-Y) IN THE TNTFRVAL ( A , B) 

XI- A 
X 2= B 

CALL FTJNCT (XI ,FX1 , DFX, INF) 

10 DO 30 1=1,20 
X= (X1 + X2J/2. 

CALL FUNCT (X,FX,DFX,INF) 

IF ( (FX1-Y) * (FX-Y) . GT. 0.) GO TO 20 
X2= X 
GO TO 30 
20 X1= X 
FX 1= FX 
30 CONTINUE 

IF (ABS (Y-FX) .LT.TOLERY) RETURN 

WRITE (6, 1000) A , B , Y 

STOP 

1000 FORMAT (////4X,49HROOT HAS FAIIEE TO CONVERGE IN THE GIVEN TNTFRVAL 
1/4 X, 3HA =,G14.6, 10X,3HB = , G 1 4 . 6 , 1 OX , 3H Y =,G14.6) 

END 
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$ I E FTC INTG1 DECK 

PEAL FUNCTION INTI (XX) 


CO PPON/GAPPP/GPP1 , GPP1 

CO PPON/C 1/GAP, R , PTZ , TTZ , UPM A CH , NST,NVP, NTURB,KPVFj,KEM,KSPTH, 

1 KS^LN , KLE,K ATCH , CTHFT , DLAP,TLAP, P^URO ,TTUPB,KPPE,KGFAD, KSDE , KLAP , 
2KPAIN,KPROF,X(100) , Y (100) ,PPFS(100) ,UE(100) , PE (100) ,POPTZ(100) , 

3 VO VCR (100) , T W A L (100) , E m A (100) 

CON MON/C 3/XCP (1O0),YOM(1OQ),S(1OO),SOL(1OO) , AF (100) ,TSE ( 100) , 

IT AWL (100) ,TAWT (100) ,TBAF (100) , PW (100) ,SW (100) ,SUTHL (100) , 

2FHSW (100) ,RHSE (100) , PE A DM (1 00) , HEA DE (100) ,NUK (100) ,KUBAP(10 0) , 

3AA (100) , PB (100) , FP (100) , DUOS (100) ,DKDS (100) , DMDL (100) 

CO.NMON/C5/SRAPL (100) ,SHAPK(100) ,B,NS 


kEAL PE,N!JW,PUBAP,tNT1 

CALL LGRNGE (SOL, ME, N ST, XX, A NS) 

INTI = ANS** (B-1 .)/ ( (1- +GMM1/2. *ANS**?) ** 

1 ( (3.+GAK— 1 . ) / (2. *GPM1) ) ) 

RETURN 

END 

$ I E FT C I NTG 2 DECK 

REAL FUNCTION INT2(XX) 

COP PON/Cl /GAP, R, PTZ, TTZ, WPP AC H, NST,NVP, NTUPB,KPVP, KFM, KSMTH, 
1KSPLN, KLE, KATCH, CTH RT,DLAP r TLAP,DTURB,7TUFB,KPRE,KGFAD,KSDE,KLAP, 
2KPAIN,K PROF, X(100),Y (IOC) ,PPES(100) , UE (100) , PE (100) , PQDTZ(IOO) , 
3VOVCR (100) ,TWAL (100) , ETA (100) 

CO P PON/C 3/XCM (100) ,YOP (100) ,S (10O) , SOL (100) , AE (100) ,TS2 (100) , 

IT AWI, (100) ,T AWT (100) , TEAR (100) ,RW (100) ,SW (100) ,SUTHL(100) , 

2RHSW (100) , RHSE (100) ,HEADW (100) ,HEABE (100) ,N0W (100) , HUBAR (100) , 

3 A A (100) , PB (100) , FF (100) ,DUDS(100),DPDS (100) ,DPDL(100) 

COPMON/C5/SB A PL (100),SHAPK(100) , B,NS 
PEAL ME,NUW,PU8AR,INT2 

IF (NS.LT.4) GO TO 10 

CALL LGRNGE (SOL, SH A FK, N S , X X , INT 2) 

RETURN 

10 DO 20 J = 2 , NS 

IF (SOL ( J ) . IT. XX) GO TO 20 

I NT 2= SHAPK ( J- 1 ) + (SHARK (J)— SHARK (J-1 ) ) * { X X- SOL ( J- 1 ) ) / (S OL ( J) -SOL (J 
1 - 1 ) ) 

RETURN 
20 CONTINUE 
RETURN 
END 


Lewis Research Center, 

National Aeronautics and Space Administration, 
Cleveland, Ohio, September 8, 1971, 
113-34. 
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APPENDIX A 


PROGRAM CHANGES FOR A GAS OTHER THAN AIR 


The program gas properties are set up for air. This program can be easily changed 
so that it applies to gases other than air. The changes are all made in subroutine 
PRECAL. 

The coefficients read in by the DATA statements of CMU, CPR, and CTC arrays 
must be changed. The equations for the curve fits have the following form: 


^sl 


a l 




Pr = a 2 





4 


k 


k 


si 


= a € 





4 


If the number of coefficients changes from five in any case, this must be reflected both 
in the DATA statements and later in the calls on CURVFT where Clvru, CPR, and CTC 
are used. If the properties are put in a different form, these cards must be removed. 

The sea-level reference values in U. S. customary and SI units of temperature 
(TSLE, TSLM), viscosity (MUSLE, MUSLM), and thermal conductivity (TCSLE, TCSLM) 
must be changed. 

The value of Sutherland's constant (TCON) and the computation of k (SUTHL (1)) 
will have to be changed. A temperature-viscosity law of the following form was used: 


JL = 

^o 



Where k gu for air was used: 

k ... / T y/Y T o + tcon \ 

SU \ T o/ V T + TCON / 

Where TCON is Sutherland's constant for air. 
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APPENDIX B 


ADDITIONAL OUTPUT 

This output is obtained when KPRE, KGRAD, KSDE, KLAM, and KPROF are set 
equal to 1. 

The output corresponding to KPRE are the geometric variables: 

X X-coordinate, m; ft 

Y Y- coordinate, m; ft 

S surface length, x, m; ft 

XOM ratio of X to C* 

YOM ratio of Y to C* 

SOL ratio of surface length to total arc length 

The next part of the output gives the local speed of sound and several temperatures: 

AE local free- stream speed of sound, m/sec; ft/sec 

TSE static temperature, K; °R 

TWAL wall temperature, K; °R 

TAWL laminar recovery temperature, K; °R 

TAWT turbulent recovery temperature, K; °R 

TBAR reference temperature, K; °R 

The final part of this output gives 

RW Reynolds number at wall, RW = (UE)(S)/NUW 

SW temperature function at wall 

SUTHL value of the coefficient in Sutherland's viscosity temperature 

o o 

RHSW static density based on the wall temperature, kg/m ; slug/ft 

3 3 

RHSE static density based on free- stream temperature, kg/m ; slug/ft 

HEADW velocity head based on the density at wall, N/m^; lbf/ft^ 

HEADE velocity head based on the free- stream density, N/m ; Ibf/ft 
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Q 2 

NUW kinematic viscosity at wall, m /sec; ft /sec 

2 

MUBAR dynamic viscosity based on reference temperature, (N)(sec)/m ; 

(lbf)(sec)/ft 2 

Output corresponding to KGRAD contains the three gradients of velocity and Mach 
number along the surface computed by finite difference methods: 

DUDS dUE/dx, sec " 1 

DMDS dME/dx, m" 1 ; ft ” 1 

DMDL dME/dSOL 

Output corresponding to KSDE contains the numerical solution of the laminar and 
turbulent differential boundary -layer equations. In the laminar case, the solution is the 
correlation number CORLN. In the turbulent case, the solution is the incompressible 
form factor FORME, and a function F of the momentum thickness. These solutions are 
printed with respect to the surface length S. 

Output corresponding to KLAM contains the variables used in the laminar subroutine 
to check for the position of instability and transition. The three variables, RTHI (in- 
creasing from station to station) and RCRIT and RTRAN (decreasing from station to sta- 
tion), are used in this analysis. When RTHI becomes larger than RCRIT, instability has 
occurred. When RTHI becomes larger than RTRAN, transition is assumed to occur. 

The variables listed are 


CORLN 

SHEAR 

DTH 

FORMTR 

SHAPL 

RTHI 

SHAPE 

RCRIT 

KBAR 

DIFF 

RTRAN 


correlation number 
shear parameter 

ratio of transformed displacement thickness to transformed momentum 
thickness 

transformed form factor 

Pohlhausen shape factor based on boundary-layer thickness 
incompressible momentum -thickness Reynolds number 
dimensionless shape factor based on momentum thickness 
critical incompressible momentum -thickness Reynolds number 
mean shape factor based on momentum thickness 

difference between transition and instability momentum -thickness Reynolds 
numbers 

incompressible momentum-thickness Reynolds number used in checking for 
transition point 
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Output corresponding to KPROF contains the velocity profiles at each station along 
the surface. The output listed is 


Y/DELTA 

Y 

Y/XMAX 

U 

U/UE 


ratio of distance normal to surface in y-direction in boundary- layer profile 
to boundary- layer thickness 

distance normal to surface in y-direction in boundary-layer profile, m; ft 
ratio of Y to XMAX 

velocity within boundary layer, m/sec; ft/ sec 
ratio of U of free- stream velocity 
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